# High-speed generation of vector beams through random spatial multiplexing

Xiao-Bo Hu1, Si-Yuan Ma1, Carmelo Rosales-Guzmán1,2 1 Wang Da-Heng
Collaborative Innovation Center, Heilongjiang Provincial Key Laboratory of
Quantum Manipulation and Control, Harbin University of Science and Technology,
Harbin 150080, China.  2 Centro de Investigaciones en Óptica, A.C., Loma del
Bosque 115, Colonia Lomas del campestre, C.P. 37150 León, Guanajuato, Mexico
carmelorosalesg<EMAIL_ADDRESS>
###### Abstract

Complex vector modes have become topical of late due to their fascinating
properties and the many applications they have found across a broad variety of
research fields. Even though such modes can be generated in a wide variety of
ways, digital holography stands out as one of the most flexible and versatile.
Along this line, Digital Micromirror Devices (DMDs) have gain popularity in
recent time due to their high refresh rates, which allows the generation of
vector modes at kHz rates. Nonetheless, most techniques are limited either by
the diversity of vector modes that can be generated or by the speed at which
they can be switched. Here we propose a technique based on the concept of
random encoding, which allows the generation of arbitrary vector beams at
speeds limited only by the refresh rate of the DMD. Our technique will be of
great relevance in research fields such as optical communications, laser
material processing and optical manipulation, amongst others.

Keywords:Complex vector beams, Digital Micromirror Devices

##  1 Introduction

Complex vector light fields, also known as vector or classically-entangled
modes, have captured the attention of many researchers across a vast variety
of research fields, such as optical tweezers, super resolution microscopy,
classical and quantum communications, optical metrology, amongst others [1, 2,
3, 4, 5, 6]. One of the main features of vector modes is their non-homogeneous
polarisation distributions, which arises from the non-separability between the
spatial and polarisation degrees of freedom [7, 8, 9, 10]. Importantly, their
mathematical structure is identical to that of quantum entangled states, one
of the main reasons why the term classically-entangled modes has been adopted
to refer to these modes [11, 12, 13]. This mathematical similarity allows to
establish many connections between the classical and the quantum world,
connections that have been the subject of an extensive research that includes
a wide variety of applications [14, 15, 16, 17, 18, 19, 20, 21]. This has, in
part, boosted the need for novel techniques capable to generate classically-
entangled modes in easier ways, with more flexibility or with high switch
speeds. Even though there are already a wide variety of techniques (see for
example [1, 22]), not all of them offer the flexibility and speed that is
required in applications such as optical communications or laser material
processing. This is the case of Spatial Light Modulators (SLMs), which offer a
great flexibility for the generation of arbitrary vector modes, but are
limited to low refresh rates ($\approx$ 60 Hz) [23, 24, 25, 26, 27, 28, 29].
In contrast, Digital Micromirror Devices (DMDs) can also generate arbitrary
vector modes at higher refresh rates [30, 31, 32, 33, 34]. In addition, DMDs
are not sensitive to polarisation, a property that has been fully exploited
not only in the generation but also in the characterisation of vector modes
[35, 36, 37, 38, 39].

While many DMD-based techniques have been developed to generate light modes
with homogeneous polarisation distribution (scalar modes), only a few have
been proposed to generate vector modes. In one of the first approaches, which
relies on the combination of a DMD and a q-plate, a scalar mode was generated
from the DMD and subsequently passed through a q-plate to convert it into a
vector mode [33]. Alternative techniques rely on the independent manipulation
of both degrees of freedom. For this, two orthogonal spatial modes with
orthogonal polarisation components are generated along independent paths and
subsequently recombined along a common propagation path, using for example
beam displacers [31]. It is worth mentioning that this approach is limited to
the generation of Laguerre-Gaussian ($LG_{p}^{\ell}$) modes since they can be
obtained as a linear superposition of Hermite-Gaussian ($HG_{nm}$) modes[40,
41]. For example, to generate a radially polarised $LG_{0}^{1}$ vector mode, a
horizontally polarised $HG_{10}$ mode is generated along one path and a
vertically polarised $HG_{01}$ mode on the other path. In a more recent
approach, which relies on the polarisation-independence property of DMDs, the
interferometric recombination was avoided by illuminating the DMD with two
orthogonally-polarised modes impinging at different angles. Nonetheless, this
technique is also limited to the set of $LG_{p}^{\ell}$ vector modes [36, 37].
Finally, it is worth mentioning that some DMDs can also display gray-scale
holograms, at the cost of reduced refresh rates, allowing the generation of
vector modes with other symmetries, this by encoding multiplexed gray-scale
hologram [34].

Even though, research fields such as optical tweezers, free-space optical
communications, laser material processing, amongst others, will certainly
benefit from vector modes with arbitrary shapes and high speeds, current
techniques does not meet both requirements. As such, in this manuscript, we
propose a novel technique capable to generate arbitrary vector modes with
almost any spatial shape reaching speeds limited only by the refresh rate of
the specific DMD (9.5 KHz for the DLP LightCrafter evaluation module). Our
technique relies on the concept of random spatial multiplexing, which was
proposed almost three decades ago in the context of phase-only optical filters
for optical pattern recognition [42] and recently in the context of beam
shaping with phase-only diffractive holograms [43, 44].

The manuscript is organised as follows, we will start by briefly explaining
the main idea behind the generation of scalar beams with DMDs. Afterwards, we
will focus on the concept of random spatial multiplexing and its application
in the generation of two independent scalar beams, which is key to generate
arbitrary vector beams at high speeds. Finally, we will demonstrate the
experimental generation of vector modes with elliptical coordinates.

##  2 Theoretical background

###  2.1 Beam shaping with digital micromirror devices

A DMD consist of an array of millions of micron-sized mirrors ($\approx 8\mu
m$), each of which can be turned to an OFF or ON state by tilting it
$-12^{\circ}$ or $+12^{\circ}$, respectively. In this way, when the DMD is
properly aligned, each mirror in the ON state will reflect light in the
desired direction [32, 45]. Since each mirror can only be in two states, it is
important to have in mind that one of the main requirements is to address the
DMD with binary holograms. Nonetheless, some devices can also deal with grays-
scale holograms, at the expense of reduced refresh rates by allowing the
micromirrors to oscillate between the ON and OFF states.

To start with, let us remind that binary-amplitude holograms are capable to
shape both the amplitude and phase of a complex light field
$U(x,y)=A(x,y)\exp(i\phi)$, where $\phi$ is the azimuthal angle given as
$\arctan(y/x)$. To this end, a binary amplitude hologram is appropriately
designed to carry in its first diffraction order the amplitude and phase
information. Further, the angle of the first diffraction order, and therefore
its position in the observation plane, can be controlled by adjusting the
spatial frequency of the grating. The desired amplitude $A(x,y)$ at a given
position $(x,y)$ in the observation plane is then acquired by locally varying
the width of each diffraction grating. In a similar way, the desired phase
$\exp(i\phi)$ is achieved by locally varying the lateral position of each
grating. More explicitly, the binary transmittance function $T(x,y)$ is given
by [45, 31],

$T(x,y)=\frac{1}{2}+\frac{1}{2}\text{sgn}\left\\{\cos\left[p(x,y)\right]+\cos\left[q(x,y)\right]\right\\},$
(1)

with,

$\begin{split}p(x,y)&=\phi(x,y)+2\pi(\nu x+\eta y),\\\
q(x,y)&=\arcsin\left({A(x,y)}/{A_{max}}\right).\end{split}$ (2)

Here, sgn$\\{\cdot\\}$ is the sign function and ${A_{max}}$ represents the
maximum value of $A(x,y)$. The term $2\pi(\nu x+\eta y)$ is a linear phase
ramp of spatial frequencies $\nu$ and $\eta$, that determines the angle of the
first diffraction order. Figure 1 (a) shows an example of a typical binary
hologram obtained from Eq. 1 for the specific case of Laguerre-Gaussian
($LG_{p}^{\ell}$) modes, as defined in [46], with radial index $p=1$ and
azimuthal index or topological charge $\ell=-1$. In this and all the following
figures showing a binary hologram, were inverted the color in such away that
the black color represents the mirrors in the ON state, which are the ones
that redirect light in the desired direction. When such transmission grating
is displayed on a DMD, we have to use the negative of such image. Notice how
the width of each grating (the black apertures) varies across the transverse
plane, resembling the amplitude of the $LG_{1}^{-1}$ mode. In addition, the
grating clearly shows the typical fork hologram, which is commonly used to
generate light beams with an azimuthal phase variation of the form
$\exp(-i\ell\phi)$, in this case $\ell=-1$. Figure 1 (b) schematically
illustrates, through a numerical simulation, the effect of such holograms when
illuminated with an expanded light beam. As expected, several diffraction
orders appear, five of which are shown here. As can be seen, the first
diffraction order carries the desired light field $U(x,y)=A(x,y)\exp(i\phi)$,
which for the sake of clarity is shown in Fig. 1 (c). In the following section
we will explain the concept of random spatial multiplexing applied to the
simultaneous generation of two spatial modes with a unique properties, such
as, phase, amplitude or diffraction angle.

Figure 1: (a) Example of a binary hologram for the case of an $LG_{1}^{-1}$
mode, (b) diffraction patter produced when such hologram is illuminated by a
flat wave front, here we only show the five diffraction orders. (c) Intensity
profile of the positive first diffraction order, which corresponds to the
encoded mode $LG_{1}^{-1}$.

###  2.2 Random binary encoding

As mentioned earlier, the concept of spatial random multiplexing has been used
before in the context of light beam shaping with phase-only holograms. In
particular, in [43], it was used to generate two Airy-vortex beams
simultaneously. Here, we will implement a similar approach where the concept
of spatial random multiplexing will be applied to multiplex two spatial modes
on a single binary amplitude hologram. To start with, we define a random
binary mask $R(x,y;a)$ having the same resolution as the DMD (1920$\times$1080
pixels), which takes values 0 or 1 at random pixels, as schematically shown in
the left panel of Fig. 2(a), for a 10$\times$10 mask. Here, $a\in[0,1]$ is a
parameter that controls the amount of pixels in the ON state. For example,
when $a=1$ all the pixels are in the ON state, whereas for $a=0.5$ only half
of the them are in the ON state. In the example shown in Fig. 2(a) only one
third of the pixels are in the ON state ($a=0.32$). When such random mask is
displayed on the DMD, it has the effect of tilting the micromirrors to the ON
or to the OFF state. Remember that only those in the ON state will diffract
light in the desired direction (the pixels represented by black color). In
this way, when a transmittance function $T(x,y)$ (see middle panel of Fig. 2
(a)) is convoluted with the random binary mask $R(x,y;a)$, some of its pixels
in the ON state are randomly switched to the OFF state (see right panel of
Fig. 2 (a)). To further exemplify this, let us take the specific case of
$LG_{p}^{\ell}$ modes, which is illustrated in Fig. 2 (b). The top panels show
the effect of multiplying the transmittance function associated to an
$LG_{2}^{2}$ mode by a random binary mask $R(x,y;a)$, as function of $a$, here
we show the cases $a$=1, 0.8, 0.6 and 0.2. Notice that a decrease of the
parameter $a$ is associated to a decrease in the number of micromirrors that
contribute to the formation of the beam in the first diffraction order. This
decrease results in a decrease of its power, as illustrated in the bottom
panels of Fig. 2 (b). Hence the parameter $a$ provides a way to directly
control the power of the mode associated to the transmittance function
$T(x,y)$. In the following section, we will further explore this idea in the
simultaneous generation of two light beams with different powers.

Figure 2: (a) A random binary mask (left) and a transmittance function
(middle) are convoluted to obtain the the mask shown on the right. (b)
Transmittance function (top) and associated intensity profile (bottom)
obtained for an $LG_{2}^{2}$ for the specific values $a=1$ $a=0.6$, $a=4$ and
$a=0.2$.

###  2.3 Spatial random multiplexing

We will now explain how the approach described above allows the simultaneous
generation of two spatially independent modes, each with unique amplitude,
phase, diffraction angle and power content. We will start by defining the
complementary function $\bar{R}=1-R(x,y;1-a)$ with the property
$R(x,y;a)\bar{R}(x,y;1-a)=0$. This implies that for $a=1$, all the pixels of
$R(x,y;1)$ are in the On state, whereas for $\bar{R}(x,y;0)$ all are in the
OFF state. On the contrary, for $a=0$, all the pixels of $R(x,y;0)$ are in the
OFF state, while for $\bar{R}(x,y;1)$ all are in the ON state. The next step
is to perform the convolution of $R(x,y;a)$ with $T_{1}(x,y)$ and
$\bar{R}(x,y;1-a)$ with $T_{2}(x,y)$ and multiplex this two holograms into a
single one. This mathematically expressed as,

$T_{f}(x,y)=R(x,y;a)*T_{1}(x,y)+\bar{R}(x,y;1-a)*T_{2}(x,y).$ (3)

This transmittance function will generate two modes simultaneously with some
of the pixels contributing to the generation of the mode associated to
$T_{1}(x,y)$ and the rest to the mode encoded in $T_{2}(x,y)$. It is important
to mention that even when some of the pixels of $T_{1}(x,y)$ overlap with some
of the pixels of $T_{2}(x,y)$, this overlaps disappears after the convolution
with the binary masks $R(x,y;a)$ and $\bar{R}(x,y;1-a)$. The previous
description is schematically illustrated in Fig.3 for the specific case of
$LG_{p}^{\ell}$. The first column shows an example of the transmittance
function $T_{1}(x,y)$ (top) and $T_{2}(x,y)$ (bottom). The second column show
an example of a random binary mask $R(x,y;a)$ (top) and its complementary
random mask $\bar{R}(x,y;1-a)$ (bottom). The inset shows an enlarged portion
of both masks to emphasize one is the complement of the other. The third row
shows the convolution $T_{1}(x,y)*R(x,y;a)$ (top) and
$T_{2}(x,y)*\bar{R}(x,y;1-a)$ bottom, which are finally added to produce the
transmittance function $T_{f}(x,y)$ shown on the right most inset of the same
figure. For this specific example, we encoded the $LG_{2}^{-1}$ mode in the
function $T_{1}(x,y)$ and $LG_{1}^{1}$ mode in the function $T_{2}(x,y)$.

Figure 3: Conceptual representation of spatial random encoding. The left
column illustrates an example for the transmittance function $T_{1}(x,y)$
(top) and $T_{2}(x,y)$ (bottom), encoding the modes $LG_{2}^{-1}$ and
$LG_{1}^{1}$. The second column show an example of a random binary mask
$R(x,y;a)$ (top) and its complementary random binary mask $\bar{R}(x,y;1-a)$
(bottom), the small inset show an enlarged portion of the mask to emphasize
that one is the complement of the other. The third column shows the
convolution $T_{1}(x,y)*R(x,y;a)$ (top) and $T_{2}(x,y)*\bar{R}(x,y;1-a)$
(bottom). The last panel shows the resultant transmittance function
$T_{f}(x,y)$.

We further illustrate this technique by using as an example the set of
$LG_{p}^{\ell}$ modes. More specifically, by generating simultaneously the
modes $LG_{2}^{-1}$ and $LG_{1}^{1}$. First, in Fig. 4 (a) we show the
transmittance function $T_{f}(x,y)$ as function of the parameter $a$ for
$a=1,0.8,0.5,0.2$ and $0$. In a similar way, in Fig. 4 (b) we show the
intensity distribution of the two generated modes. Notice that for $a=1$ only
the hologram corresponding to the $LG_{2}^{-1}$ mode appears, as $a$
decreases, the hologram corresponding to the mode $LG_{1}^{1}$ starts
appearing, while the first one gradually vanishes until it completely
disappears for $a=0$. In Fig. 4 (b) we show the far-field intensity profile
produced by the holograms shown in Fig. 4 (a), where the first diffraction
order is shown and the 0th is only indicated for reference purposes. As
expected, for $a=1$ only the $LG_{2}^{-1}$ mode is present, as we decrease
$a$, its power decreases while the power of the $LG_{1}^{1}$ mode increases,
for $a=0.5$ the power of both beams is the same. Finally for $a=0$, only the
mode $LG_{1}^{1}$ is present.

Figure 4: (a) Transmittance function $T_{f}(x,y)$ and (b) Far field intensity
distribution of as function of the parameter $a$ for $a=1$, 0.8, 0.5, 0.2 and
0. The examples shown correspond to the modes $LG_{2}^{-1}$ (left) and
$LG_{1}^{1}$ (right).

###  2.4 Generation of vector modes through spatial random multiplexing

To start this section, let us remind that vector modes result from a non-
separable coupling between the spatial and polarisation degrees of freedom.
Experimentally, such modes can be generated as a coaxial superposition of two
orthogonal spatial modes carrying orthogonal polarisation. Mathematically,
this can be expressed as [7, 9]

${\bf U}({\bf r})=\cos(\theta)U_{1}({\bf r})\hat{\bf
e}_{1}+\sin(\theta)\rme^{i\alpha}U_{2}({\bf r})\hat{\bf e}_{2},$ (4)

where $U_{1}({\bf r})$ and $U_{2}({\bf r})$ are two normalized complex-valued
orthogonal functions representing the spatial degree of freedom, whereas
$\hat{\bf e}_{1}$ and $\hat{\bf e}_{2}$ are two orthogonal vectors
representing the state of polarisation. The coefficient $\theta\in[0,\pi/2]$
allows a transition between scalar ($\theta=0$ and $\theta=\pi/2$) and vector
($\theta=\pi/4$). The complex term $\exp(i\alpha)$, with $\alpha\in[0,\pi]$,
represents an inter-modal phase between both constituting spatial modes.

We will explain how to use random binary encoding to generate vector beams
with arbitrary spatial shape and polarisation distribution. The main idea
behind this technique relies on the ability to generate two independent and
orthogonal spatial modes whose power content can be controlled also
independently. The requirement of having two independent modes is essential as
we also need to manipulate their polarisation state independently, as required
by Eq. 4. Once generated, these modes are recombined along a common
propagation axis where the desired vector beam is generated. To achieve this
superposition we will use a recently proposed method, which consist on
illuminating the DMD with two modes of orthogonal polarisation impinging at
different angles but exactly on the same spot, which coincides with the
geometric center of the binary hologram. On the DMD we display a multiplexed
hologram encoding two orthogonal spatial modes. Here is where our random
spatial multiplexing technique becomes crucial since it allows us to generate
two independent modes with controllable power using the transmittance function
$T_{f}$ given by Eq. 3. In this way, when the two beams with orthogonal
polarisation impinge on the DMD, two copies of the modes encoded in $T_{f}$
are generated in the first diffraction order, one with one polarisation, lets
say right circular, and the other with the orthogonal polarisation, in this
case left circular. All we have to do now is to adjust the grating periods of
each mode to overlap two of the four modes along the same propagation axis.
More details of this procedure will be provided in the following section.

##  3 Experimental generation of vector modes through spatial random
multiplexing

###  3.1 Experimental setup

To generate the modes described by Eq. 4 we will use the approach described in
[36], based on the experimental setup depicted in Fig. 5(a). Here, a
horizontally polarised laser beam ($\lambda=532$ nm, 300 mW) is first expanded
with the pair of lenses L1 and L2 of focal lengths $f_{1}=20$ mm and
$f_{2}=200$ mm, respectively. The beam is then sent through a Half-wave plate
(HWP) oriented at $22.5^{\circ}$, to obtain a diagonally polarised beam.
Afterwards it is split into its horizontal and vertical polarisation
components using a Polarising Beam Splitter (PBS). Both beams are then
directed to a polarisation-independent Digital Micromirror Device (DMD, DLP
Light Crafter 6500 from Texas Instruments), impinging at slightly different
angles ($\approx 1.5^{\circ}$) at the center of the DMD. Here, a spatial
random multiplexed hologram is displayed, which is generated as explained
before using the transmittance function given by Eq. 3. The hologram contains
the two constituting wave fields $U_{1}(x,y)$ and $U_{2}(x,y)$ required to
generate the vector mode given by Eq.4, which is shown in Fig. 5(b). As
result, four modes will appear in the first diffraction order, two for each
beam impinging on the DMD, as shown in Fig. 5(c), where modes m1 and m2
correspond to Beam 1 and modes m3 and m4 to Beam 2. Crucially, modes m2 and
m3, which are already orthogonal, also carry orthogonal polarisation, as
indicated in Fig. 5(c). Hence, all we have to do is to adjust the period of
the diffraction grating of each beam to overlap them along a common
propagation axis. The resulting mode is already the desired vector beam, which
we isolate with the help of a Spatial Filter (SF), located at the focusing
plane of a telescope composed by lenses $\rm L_{3}$ and $\rm L_{4}$, both with
focal lengths $f=100$ mm. A Quarter-Wave plate (QWP1) can be added either
before or after the DMD to change the vector mode from the linear to the
circular polarisation basis. Finally, we add a lens $\rm L_{5}$ of focal
length $f=300$ mm to send the generated modes to a Charge-Coupled Device (CCD)
camera (FLIR FL3-U3-120S3C-C from Pointgray).

Figure 5: (a) Experimental setup implemented to generate vector modes using
binary random encoding. A horizontally polarized laser beam, expanded and
collimated with lenses L1 and L2 is changed to diagonal polarisation with a
Half-Wave Plate (HWP1) oriented at 45∘ and separated afterwards into its two
polarisation components using a Polarising Beam Splitter (PBS). Both beams are
then sent to the centre of a Digital Micromirror Device (DMD), where the
hologram is displayed. A Spatial Filter (SF) placed in the focal plane of a
telescope formed by lenses L3 and L4 filters the desired order. A Quarter-Wave
Plate (QWP1) at 45∘ ensures the vector mode is in the circular polarisation
basis, which is imaged with lens L5 to a Charge-Coupled Device (CCD) camera.
The state of polarisation is reconstructed through Stokes polarimetry by means
of a HWP, a QWP and a polariser (P). (b) Typical hologram displayed on the
DMD. (c) Intensity profile of the first diffraction order generated by Beam1
and Beam2, which generates modes m1, m2 and m3, m4, respectively, with
orthogonal polarisation, represented by the orange and green circular arrows.

###  3.2 Polarisation reconstruction

One way to analyse the generated vector modes is by reconstructing their state
of polarisation through the Stokes parameters, which can be computed from a
series of intensity measurements as [47],

$\begin{split}\@add@centering\centering&S_{0}=I_{0},\hskip
54.06023ptS_{1}=2I_{H}-S_{0},\hskip 2.84526pt\\\ &S_{2}=2I_{D}-S_{0},\hskip
28.45274ptS_{3}=2I_{R}-S_{0},\end{split}$ (5)

Here $I_{0}$ is the total intensity of the beam, $I_{H}$, $I_{D}$ and $I_{R}$
are the corresponding intensities of the horizontal, diagonal and right-handed
polarisation components. Experimentally, these intensity are measured using a
series of polarisation filters, composed by a linear polariser (P), a HWP and
a QWP, as indicated in Fig. 5, and recorded with the CCD camera [39, 38]. More
precisely, the intensities of the horizontal ($I_{H}$) and diagonal ($I_{D}$)
polarisation components, are obtained by inserting a linear polariser at
$0^{\circ}$ and $45^{\circ}$, respectively, whereas the intensity
corresponding to the right circular polarisation component ($I_{R}$) is
acquired by inserting a QWP at $45^{\circ}$ in combination with a linear
polariser at $90^{\circ}$.

###  3.3 Experimental generation of Laguerre-Gaussian vector beams

As a first example we will show the arbitrary generation of cylindrical vector
modes encoded in the $LG_{0}^{\ell}$ spatial modes. For the sake of clarity,
we rewrite Eq. 4 using this spatial basis and the circular polarisation basis,
$\hat{\bf e}_{R}$ and $\hat{\bf e}_{L}$,

${\bf U}_{LG}({\bf r})=\cos\theta LG_{p_{1}}^{\ell_{1}}({\bf r})\hat{\bf
e}_{R}+\sin\theta\rme^{i\alpha}LG_{p_{1}}^{\ell_{2}}({\bf r})\hat{\bf e}_{L}.$
(6)

Examples of the typical vector beams that can be generated are shown in Fig. 6
for the specific cases $\\{LG_{0}^{1};LG_{0}^{2}\\}$ (Lemon),
$\\{LG_{0}^{2};LG_{0}^{1}\\}$ (Star), $\\{LG_{0}^{-2};LG_{0}^{2}\\}$ (Spider)
and $\\{LG_{0}^{2};LG_{0}^{-2}\\}$ (Web). An example of the multiplexed
holograms required to generate these modes is shown in Fig. 6 (a). In a
similar way, Fig. 6 (b) shows a theoretical simulation of the transverse
intensity profile, overlapped with the transverse polarisation distribution.
Finally, Fig. 6 (c) shows the experimental intensity profile, also overlapped
with the reconstructed polarisation distribution. Small deviation between the
numerical simulations and our experimental results are to be expected since
all DMDs introduce optical aberrations caused by the non-flatness of the
screen, nonetheless, they can be minimized by adding a correction mask [32].

Figure 6: Experimental generation of the Laguerre-Gaussian vector modes
$\\{LG_{0}^{1};LG_{0}^{2}\\}$ (Lemon), $\\{LG_{0}^{2};LG_{0}^{1}\\}$ (Start),
$\\{LG_{0}^{-2};LG_{0}^{2}\\}$ (Spider) and $\\{LG_{0}^{2};LG_{0}^{-2}\\}$
(Web). (a) Examples of the holograms displayed on the DMD, where, for the
purpose of display, we inverted the colors. (b) Numerical simulation of
transverse intensity profile overlapped with its corresponding polarisation
distribution. (c) Experimentally generated modes.

As a second example, we show the case of Laguerre-Gaussian vector beams with
arbitrary values of $\theta$ and $\alpha$. Without loss of generality, we
restrict this example to the vector mode formed by the superposition of the
modes $LG_{1}^{2}$ and $LG_{1}^{-2}$. As mentioned earlier, $\theta$ modulates
the contribution of each spatial mode, which we achieved experimentally by
varying the coefficient $a$ in the transmittance function associated to each
mode, namely, $R(x,y;a)*LG_{1}^{2}$ and $\bar{R}(x,y;1-a)*LG_{1}^{-2}$. The
intensity profile of modes generated for increasing values of $a\in[0,1]$
($\theta\in[0,\pi/2$]), are shown as a media file ( Media 1), selected frames
this video are shown in Fig. 7 (a). In order to see the effect produced when
$\theta$ is increased, we passed the vector beam through a linear polariser
oriented horizontally. As result of this, the intensity profile of the mode
observed with the CCD, gradually changes from two concentric rings of light
(scalar), to a petal-like structure (vector) and back to the structure of
concentric rings (scalar). Finally, the effect of changing the inter-modal
phase associated to $\alpha$ can be achieved experimentally, by simply adding
an extra phase to one of the constituting modes and encoded directly on its
transmittance function. Again, the effect of $\alpha$ can be clearly observed
by sending the vector beam through a linear polariser. As result, the
intensity distribution features a concentric structure of four petals that
will rotate anticlockwise with increasing values of $\alpha$. Figure 7 shows a
selected set of four frames extracted from  Media 2, where this effect is
clearly shown.

Figure 7: (a) Selected frames taken from  Media 1 where the intensity
distribution of the vector modes as function of $\theta\in[0,\pi/2]$ is shown.
(b) Selected frames taken from  Media 2 where the intensity distribution as
function of $\alpha\in[0,\pi]$ is shown. In both cases, the vector beam was
sent through a linear polarise oriented horizontally.

###  3.4 Experimental generation of elliptical Ince-Gaussian beams

Finally, and to show the capabilities of our technique, in this section we
demonstrate the generation of vector modes with elliptical symmetry by
encoding the spatial degree of freedom in the Helical Ince-Gaussian (HIG)
vector modes. Such vector modes have been previously generated with phase
holograms encoded in SLMs [48] and DMDs [34] and therefore with low refresh
rates. Importantly, our technique allows us to generate such modes using
binary holograms at speeds only limited by the specific characteristic of the
DMD, in our case 9.5 kHz. Let us first remind that IG modes are natural
solutions of the paraxial wave equation in elliptical cylindrical coordinates
${\bf r}=(\xi,\eta,z)$, with, $\xi\in[0,\infty)$ and $\eta\in[0,2\pi)$
representing the radial and angular elliptical coordinates, respectively [49].
Here, we will not provide specific details of how this solutions are derived,
but only the mathematical expressions that describe them, more details can be
found in [49]. The scalar IG mode can be classified into even or odd in pure
relation to the even and odd Ince polynomials, $C_{p}^{m}(\cdot,\varepsilon)$
and $S_{p}^{m}(\cdot,\varepsilon)$, which constitute these modes.
Mathematically, they have the specific form,

$\begin{split}IG_{q,m,\varepsilon}^{e}({\bf
r})&=\frac{C\omega_{0}}{\omega(z)}C_{q}^{m}(i\xi,\varepsilon)C_{q}^{m}(\eta,\varepsilon)\text{e}^{-\frac{r^{2}}{\omega(z)}}\text{e}^{-i\left(kz+Z-\Phi\right)},\\\
IG_{q,m,\varepsilon}^{o}({\bf
r})&=\frac{S\omega_{0}}{\omega(z)}S_{q}^{m}(i\xi,\varepsilon)S_{q}^{m}(\eta,\varepsilon)\text{e}^{-\frac{r^{2}}{\omega(z)}}\text{e}^{-i\left(kz+Z-\Phi\right)},\end{split}$
(7)

where $C$ and $S$ are normalisation constants and the superscripts $e$ and $o$
stand for even and odd, respectively. Further, the subscripts
$p,m\in\mathbb{N}$, which are restricted by $0\leq m\leq q$ for even functions
and $1\leq m\leq q$ for odd functions, provide the order of the mode.
Additional parameters, are related to the Gaussian envelope of beam waist
$\omega_{0}$ that confines these modes. More precisely
$\Phi=(p+1)\arctan(z/z_{R})$ is the Gouy phase, $Z(z)=kr^{2}/2R(z)$ is an
additional phase term related to the radius of curvature $R(z)=z+z_{R}^{2}/z$
of the phase front and $z_{R}=\pi\omega_{0}^{2}/\lambda$ is the Rayleigh
length. In addition, the coherent superposition of the even and odd modes
gives rise to the HIG modes $IG_{q,m,\varepsilon}^{h+}({\bf r})$ and
$IG_{q,m,\varepsilon}^{h-}({\bf r})$, which can be expressed mathematically
as,

$\begin{split}IG_{q,m,\varepsilon}^{h+}({\bf
r})&=IG_{q,m,\varepsilon}^{e}({\bf r})+iIG_{q,m,\varepsilon}^{o}({\bf r}),\\\
IG_{q,m,\varepsilon}^{h-}({\bf r})&=IG_{q,m,\varepsilon}^{e}({\bf
r})-iIG_{q,m,\varepsilon}^{o}({\bf r}),\end{split}$ (8)

Perhaps, one of the most salient properties of IG modes is that they represent
a larger family of modes that includes the Laguerre- and Hermite-Gaussian
modes as specific cases, which is related to the ellipticity parameter,
$\varepsilon=2f_{0}/\omega_{0}^{2}$, $\varepsilon\in[0,\infty)$. More
precisely, the $IG_{q,m,\varepsilon}^{h}$ transform into the Laguerre-Gaussian
($LG_{p}^{\ell}$) when $\varepsilon=0$, with their indices related as $\ell=m$
and $p=(q-m)/2$. In a similar way, the $IG_{p,m,\varepsilon}^{e,o}$ modes
become the Hermite Gaussian ($HG_{n_{x}n_{y}}$) modes for
$\varepsilon\to\infty$, with their indices related as $n_{x}=m$, $n_{y}=p-m$,
for even modes and $n_{x}=m-1$, $n_{y}=p-m+1$ for odd modes [50]. Notably, the
set of $IG_{q,m,\varepsilon}^{h+}({\bf r})$ modes are orthogonal to the set of
modes $IG_{q,m,\varepsilon}^{h-}({\bf r})$. Hence, we can use these as the
spatial degree of freedom to generate HIG vector modes. With this restriction,
Eq. 4 becomes,

${\bf U}_{IG^{h}}({\bf r})=\cos\theta IG_{q,m,\varepsilon}^{h+}({\bf
r})\hat{\bf e}_{R}+\sin\theta\rme^{i\alpha}IG_{q,m,\varepsilon}^{h-}({\bf
r})\hat{\bf e}_{L}$ (9)

To begin with, we will show that our technique allows the generation of
arbitrary IG vector modes as represented on a High Order Poincaré Sphere
(HOPS) [34]. This geometric representation, provides a way to visualize
arbitrary vector modes on the surface of a unitary sphere by associating their
inter-modal phase ($\alpha$) and weighting coefficient ($\theta$) to a point
with coordinates $(2\theta,2\alpha)$. In this way, the two scalar modes with
right and left circular polarisation are mapped to the North and South Poles,
respectively, and pure vector modes to the equator. Without loss of
generality, we restrict our analysis to the specific case given by HIG vector
modes ${\rm HIG}_{5,3,2}^{h}({\bf r})$. As a first example, Fig. 8 (a) shows a
HOPS representation of five vector modes labeled from 1 to 5, with specific
coordinates, in the same order, $(0,0),(\pi/2,0),(\pi,0),(\pi/2,3\pi/4)$ and
$(3\pi/2,\pi/4)$. In this representation, the Cartesian coordinates of the
HOPS are related to the higher-order Stokes parameters $S_{1}^{m}$,
$S_{2}^{m}$ and $S_{3}^{m}$, as detailed in [51]. The holograms required to
generate such modes, which are also computed in the same way as explained
before but with the HIG modes, are shown in the top row of Fig. 8 (a). Note
that this holograms do not correspond to the holograms displayed in the DMD,
they have been scaled for display purposes. The middle row of Fig. 8 (b) shows
a numerical simulation of the transverse intensity profile overlapped with the
polarization distribution of this specific cases. Our experimental results are
shown in the bottom row of Fig. 8 (b), showing a high agreement with our
numerical simulations.

Figure 8: (a) Geometric representation of elliptical vector modes on a high
order Poincaré sphere. (b) Typical hologram displayed on the DMD (top),
numerical simulation of the intensity and polarisation distribution (middle),
experimental results (bottom).

Our technique is also capable to generate HIG vector modes with arbitrary
values of ellipticity. As a way of example, Fig. 9 shows a representative set
of examples for the specific case ${\rm HIG}_{5,3,\varepsilon}$ and
$\varepsilon\in[0,2,5,100]$. The top row shows an example of the required
holograms to generate such modes, the middle row a numerical simulation of the
intensity and polarisation distribution and the third row our experimental
results. Notice the high similitude between both sets.

Figure 9: Helical Ince-Gaussian vector modes as a function of the ellipticity
$\varepsilon$. Binary hologram(top row) numerical simulations (middle row) and
experimental results (bottom row) of the transverse intensity profiles
overlapped with the polarization distribution. Here, we show the specific case
${\rm IG}_{5,3,\varepsilon}^{h}$ for $\varepsilon=0,2,5$ and 100. Notice the
transition from Laguerre- to Hermite-Gaussian vector modes.

##  4 Conclusions

Vector beams represent a general state of light in which the spatial and
polarisation degrees of freedom are coupled in a non-separable way. The unique
properties of these modes have become of great relevance in applications such
as optical communications, laser material processing, optical metrology,
amongst others. Nonetheless, most applications are based on vector modes with
cylindrical symmetry, such as Laguerre- or Bessel-Gaussian modes. Noteworthy,
the availability of such modes with other symmetries, such as elliptic or
parabolic, will benefit current applications and pave the path towards the
development of new. In addition, the generation of vector modes with high
refresh rates provide an additional advantage that will certainly enhance the
capabilities of these. Hence, it is almost natural to seek for novel technique
which can offer high flexibility and high generation rates. Nonetheless,
existing techniques capable to generate vector beams with arbitrary spatial
shape are limited to low refresh rates and those with high generation rates
are limited to the reduced set of cylindrical vector modes. Here, we proposed
a technique that offers both advantages simultaneously. This technique is
based on the concept of random binary spatial multiplexing, which in
combination with the Digital Micromirror Device (DMD) technology, is capable
to generate vector modes with arbitrary spatial distribution, at the maximum
generation speed offered by the specific DMD model. As we demonstrated, random
spatial multiplexing allow us to generate two scalar modes with independent
properties, such as spatial shape, phase profile, polarisation, power content
and diffraction angle. The main idea consist on using some pixels of the DMD,
selected at random, to generate one scalar mode and the remaining ones to
generate the other mode. These two modes are then recombined along a common
propagation axis, using a polarisation-insensitive approach [36], where the
vector mode is generated. Our approach provides with a powerful tool for the
generation of vector modes with arbitrary spatial shapes and high refresh
rates, which will be of great relevance in various applications. For example,
in the field of optical tweezers, for the localized tailoring of optical
forces [4], or in free-space optical communications for the generation of
vector modes with high resilience to turbulence [22] and last but not least,
in laser material processing to enhance the drilling capabilities of the laser
beams employed [52].

## Acknowledgement

The authors would like to thank Dr. Valeria Rodriguez-Fajardo for fruitful
discussions.

## Funding

This work was partially supported by the National Natural Science Foundation
of China (NSFC) under Grant No. 61975047.

## Disclosures

The authors declare that there are no conflicts of interest related to this
article

## References

## References

  * [1] Rosales-Guzmán C, Ndagano B and Forbes A 2018 J. Opt. 20 123001
  * [2] Rubinsztein-Dunlop H, Forbes A, Berry M V, Dennis M R, Andrews D L, Mansuripur M, Denz C, Alpmann C, Banzer P and Bauer T 2017 J. Opt. 19 013001
  * [3] Zhan Q 2009 Adv. Opt. Photonics 1 1–57
  * [4] Bhebhe N, Williams P A C, Rosales-Guzmán C, Rodriguez-Fajardo V and Forbes A 2018 Sci. Rep. 8 17387
  * [5] Berg-Johansen S, Töppel F, Stiller B, Banzer P, Ornigotti M, Giacobino E, Leuchs G, Aiello A and Marquardt C 2015 Optica 2 864–868
  * [6] Hu X B, Zhao B, Zhu Z H, Gao W and Rosales-Guzmán C 2019 Opt. Lett. 44 3070–3073
  * [7] Galvez E J, Khadka S, Schubert W H and Nomoto S 2012 Appl. Opt. 51 2925–2934
  * [8] Beckley A M, Brown T G and Alonso M A 2010 Opt. Express 18 10777–10785
  * [9] Galvez E J 2015 Light Beams with Spatially Variable Polarization (Wiley-Blackwell) chap 3, pp 61–76 ISBN 9781119009719
  * [10] Otte E, Alpmann C and Denz C 2016 J. Opt. 18 074012
  * [11] Konrad T and Forbes A 2019 Contemporary Physics 1–22
  * [12] Forbes A, Aiello A and Ndagano B 2019 Classically entangled light Progress in Optics (Elsevier Ltd.) pp 99–153
  * [13] Toninelli E, Ndagano B, Vallés A, Sephton B, Nape I, Ambrosio A, Capasso F, Padgett M J and Forbes A 2019 Advances in Optics and Photonics 11 67–134
  * [14] Balthazar W F, Souza C E R, Caetano D P, ao E F G, Huguenin J A O and Khoury A Z 2016 Opt. Lett. 41 5797–5800
  * [15] Diego G S, Robert B, Felix Z, Christian V, Markus G, Matthias H, Stefan N, Michael D, Andrea A, Marco O and Alexander S 2016 Laser Photonics Rev. 10 317–321
  * [16] Eberly J H, Qian X F, Qasimi A A, Ali H, Alonso M A, Gutiérrez-Cuevas R, Little B J, Howell J C, Malhotra T and Vamivakas A N 2016 Phys. Scr 91 063003 URL http://stacks.iop.org/1402-4896/91/i=6/a=063003
  * [17] Li P, Wang B and Zhang X 2016 Opt. Express 24 15143
  * [18] Qian X F, Vamivakas A N and Eberly J H 2017 Opt. Photon. News 28 34–41
  * [19] Borges C V S, Hor-Meyll M, Huguenin J A O and Khoury A Z 2010 Phys. Rev. A 82(3) 033833
  * [20] Ndagano B, Perez-Garcia B, Roux F S, McLaren M, Rosales-Guzmán C, Zhang Y, Mouane O, Hernandez-Aranda R I, Konrad T and Forbes A 2017 Nature Phys. 13 397–402
  * [21] Töppel F, Aiello A, Marquardt C, Giacobino E and Leuchs G 2014 New J. Phys. 16 073019
  * [22] Ndagano B, Nape I, Cox M A, Rosales-Guzmán C and Forbes A 2018 J. Light. Technol. 36 292–301
  * [23] Dudley A, Li Y, Mhlanga T, Escuti M and Forbes A 2013 Opt. Lett. 38 3429–3432
  * [24] Chen S, Zhou X, Liu Y, Ling X, Luo H and Wen S 2014 Opt. Lett. 39 5274–5276 URL http://ol.osa.org/abstract.cfm?URI=ol-39-18-5274
  * [25] Wang X l, Ding J, Ni W j, Guo C s and Wang H t 2007 Opt. Lett. 32 3549–3551
  * [26] Maurer C, Jesacher A, Fürhapter S, Bernet S and Ritsch-Marte M 2007 New J. Phys. 9 78
  * [27] Moreno I, Davis J A, Hernandez T M, Cottrell D M and Sand D 2012 Opt. Express 20 364–376
  * [28] Rosales-Guzmán C and Forbes A 2017 How to shape light with spatial light modulators SPIE.SPOTLIGHT (SPIE Press)
  * [29] Rosales-Guzmán C, Bhebhe N and Forbes A 2017 Opt. Express 25 25697–25706
  * [30] Ren Y X, Lu R D and Gong L 2015 Annalen der Physik 527 447–470
  * [31] Mitchell K J, Turtaev S, Padgett M J, Čižmár T and Phillips D B 2016 Opt. Express 24 29269–29282
  * [32] Scholes S, Kara R, Pinnell J, Rodríguez-Fajardo V and Forbes A 2019 Optical Engineering 59 1 – 12
  * [33] Gong L, Ren Y, Liu W, Wang M, Zhong M, Wang Z and Li Y 2014 J. Appl. Phys. 116 183105
  * [34] Yao-Li, Hu X B, Perez-Garcia B, Bo-Zhao, Gao W, Zhu Z H and Rosales-Guzmán C 2020 Applied Physics Letters 116 221105
  * [35] Zhao B, Hu X B, Rodríguez-Fajardo V, Forbes A, Gao W, Zhu Z H and Rosales-Guzmán C 2020 Applied Physics Letters 116 091101 (Preprint https://doi.org/10.1063/1.5142163) URL https://doi.org/10.1063/1.5142163
  * [36] Rosales-Guzmán C, Hu X B, Selyem A, Moreno-Acosta P, Franke-Arnold S, Ramos-Garcia R and Forbes A 2020 Scientific Reports 10 10434
  * [37] Selyem A, Rosales-Guzmán C, Croke S, Forbes A and Franke-Arnold S 2019 Phys. Rev. A 100(6) 063842
  * [38] Manthalkar A, Nape I, Bordbar N T, Rosales-Guzmán C, Bhattacharya S, Forbes A and Dudley A 2020 Opt. Lett. 45 2319–2322
  * [39] Zhao B, Hu X B, Rodríguez-Fajardo V, Zhu Z H, Gao W, Forbes A and Rosales-Guzmán C 2019 Opt. Express 27 31087–31093
  * [40] Beijersbergen M, Allen L, van der Veen H and Woerdman J 1993 Optics Communications 96 123–132
  * [41] O’Neil A T and Courtial J 2000 Optics Communications 181 35–45
  * [42] Davis J A and Cottrell D M 1994 Opt. Lett. 19 496–498
  * [43] Rosales-Guzmán C, Mazilu M, Baumgartl J, Rodríguez-Fajardo V, Ramos-García R and Dholakia K 2013 Journal of Optics 15 044001
  * [44] Fuentes J L M and Moreno I 2018 Opt. Express 26 5875–5893
  * [45] Mirhosseini M, Magana-Loaiza O S, Chen C, Rodenburg B, Malik M and Boyd R 2013 Opt. Express 21
  * [46] Hu X B, Dong M X, Zhu Z H, Gao W and Rosales-Guzmán C 2020 Scientific Reports 10 199
  * [47] Goldstein D H 2011 Polarized light (CRC Press)
  * [48] Otte E and Denz C 2018 Opt. Lett. 43 5821–5824
  * [49] Bandres M A and Gutiérrez-Vega J C 2004 Opt. Lett. 29 144–146
  * [50] Bandres M A and Gutiérrez-Vega J C 2004 J. Opt. Soc. Am. A 21 873–880 URL http://josaa.osa.org/abstract.cfm?URI=josaa-21-5-873
  * [51] Milione G, Sztul H I, Nolan D A and Alfano R R 2011 Phys. Rev. Lett. 107(5) 053601
  * [52] Meier M, Romano V and Feurer T 2007 Appl. Phys. A 86 329–334


[a]P. Dimopoulos

#  $K$\- and $D_{(s)}$-meson leptonic decay constants with physical light,
strange and charm quarks by ETMC

R. Frezzotti  M. Garofalo  S. Simula

###### Abstract

We present a lattice QCD computation and preliminary results for the leptonic
decay constants of the pseudoscalar mesons $K$, $D$ and $D_{s}$ in the
isosymmetric QCD limit. The computation is based on simulations of
$N_{f}=2+1+1$ dynamical quarks performed by the Extended Twisted Mass
Collaboration (ETMC), where the light, strange and charm quark masses are all
tuned at their physical values. We also present preliminary unitarity checks
for the first and second rows of the Cabibbo-Kobayashi-Maskawa matrix.

##  1 Introduction

Of the most important hadronic inputs for obtaining estimates of the Cabibbo-
Kobayashi-Maskawa (CKM) matrix elements (ME) are the values of the leptonic
decay constants of pseudoscalar (PS) mesons. In these proceedings we present a
high precision lattice QCD (LQCD) calculation for the $K$, $D$ and $D_{s}$ PS-
meson decay constants using $N_{f}=2+1+1$ gauge ensembles generated by the
Extended Twisted Mass Collaboration (ETMC).

##  2 Lattice action and simulations

We employ the twisted mass (tm) fermionic formulation that ensures automatic
O$(a)$-improvement for all observables as far as tuning at maximal twist is in
place Refs [1, 2]. Moreover as it has been shown in Refs [3, 4, 5, 6] the
inclusion of the clover term in the maximally twisted fermionic action
provides the beneficial property of reduced O$(a^{2})$ cutoff and isospin
breaking effects.

ETMC has performed $N_{f}=2+1+1$ dynamical quark simulations employing the sea
quark action written as $S_{sea}=S_{g}^{\text{Iwa}}+S_{tm}^{\ell}+S_{tm}^{h}$,
where $S_{g}^{\text{Iwa}}$ is the Iwasaki improved gauge action [7] and

$\displaystyle S_{tm}^{\ell}$
$\displaystyle=\sum_{x}\bar{\chi}_{\ell}(x)\left[D_{W}(U)+\frac{i}{4}c_{\mathrm{SW}}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}(U)+m_{0\ell}+i\mu_{\ell}\tau^{3}\gamma^{5}\right]\chi_{\ell}(x)\,,$
(1) $\displaystyle S_{tm}^{h}$
$\displaystyle=\sum_{x}\bar{\chi}_{h}(x)\left[D_{W}(U)+\frac{i}{4}c_{\mathrm{SW}}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}(U)+m_{0h}-\mu_{\delta}\tau_{1}+i\mu_{\sigma}\tau^{3}\gamma^{5}\right]\chi_{h}(x)\,,$
(2)

are, respectively, the fermionic actions used in the light ($S_{tm}^{\ell}$)
and strange-charm ($S_{tm}^{h}$) quark sectors. Notice that in Eq. (1)
$\chi_{\ell}=(u,d)^{T}$ represents the light quark doublet, while the
degenerate light twisted and the (untwisted) Wilson quark masses are denoted
by $\mu_{\ell}$ and $m_{0\ell}$, respectively. For the heavy quark action, Eq.
(2), the doublet $\chi_{h}=(s,c)^{T}$ represents the mass non-degenerate
strange and charm quarks. Here $m_{0h}$ denotes the (untwisted) Wilson quark
mass while the parameters $\mu_{\delta}$ and $\mu_{\sigma}$ in combination
with the presence of the Pauli matrices $\tau_{1}$ and $\tau_{3}$ lead to
quark mass non-degeneracy, see Ref. [2]. In both equations $D_{W}(U)$ is the
massless Wilson-Dirac operator, and the Sheikoleslami-Wohlert improvement
term, $c_{\mathrm{SW}}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}(U)$, has been
included for the reason already explained above. The value for the clover
parameter $c_{\mathrm{SW}}$ is set by using the 1–loop tadpole boosted
estimate as presented in Ref. [8]. The condition of maximal twist is achieved
by tuning the hopping parameter $\kappa$ for the untwisted Wilson quark mass
such as $m_{0\ell}=m_{0h}=m_{crit}$. Details about the lattice action and the
algorithmic setup are presented in Refs [3, 9, 10].

Simulations have been carried out reaching the physical mass values of both
light and heavy (strange and charm) quarks. As for the latter the sea quark
mass parameters ($\mu_{\sigma}$ and $\mu_{\delta}$) have been tuned so that
the two phenomenological conditions $m_{c}/m_{s}=11.8$ and
$m_{D_{s}}/f_{D_{s}}=7.9$ [11] are accurately reproduced by each of the
$N_{f}=2+1+1$ ensembles.

In order to avoid O$(a^{2})$ mixing effects in the physical observables
involving the heavy quarks (strange and charm), owing to the form of the sea
quark action of Eq. (2), we opted for a non-unitary lattice setup. Therefore
in the valence sector we employ the Osterwalder-Seiler fermionic
regularisation [12] which treats the strange and charm quarks in a flavour
diagonal way. The valence action in the strange and charm sectors is given by:

$S_{val}^{f}=\sum_{x}\bar{\chi}_{f}^{val}(x)\left(D_{W}^{cr}(U)+\frac{i}{4}c_{\mathrm{SW}}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}(U)+i\gamma_{5}\mu_{f}\right)\chi_{f}^{val}(x)\,\,,f=s,c,$
(3)

where $D_{W}^{cr}(U)\equiv D_{W}(U)|_{m_{0}=m_{crit}}$ is the critical Wilson-
Dirac operator and $\chi_{f}^{val}$ denotes single quark flavour field. It has
been shown in Ref. [13] that this kind of mixed action preserves the automatic
O$(a)$-improvement of physical observables, i.e. lattice artifacts, including
those violating unitarity, scale as O$(a^{2})$ implying that unitarity is
safely recovered in the continuum limit.

In this study we use $N_{f}=2+1+1$ simulations performed at three values of
the lattice spacing in the range $[0.69,0.95]$ fm and at several pseudoscalar
mass values spanning from the physical pion mass up to 350 MeV. In Table 1
essential simulation details are presented. For the computation of $w_{0}/a$
and the gradient-flow $w_{0}$-determination, we refer the reader to Refs [14,
9, 15]. The scale setting is performed using the isosymmetric QCD value
$f_{\pi}^{isoQCD}=130.4(2)$ [16].

The present lattice computation is performed in the isoymmetric QCD limit.
However future work based on the same gauge ensembles is planned in order to
take into account isospin breaking effects along the lines of the work of e.g.
Ref [17].

$\beta$ | Ens. | $(L,T)$ |  $M_{\pi}$ (MeV) | # meas. | $w_{0}/a$
---|---|---|---|---|---
1.726 | cA211.12.48 | (48,96) | 167 | 322 | 1.8355(35)
cA211.30.32 | (32,64) | 261 | 1237
cA211.40.24 | (24,48) | 302 | 662
cA211.53.24 | (24,48) | 346 | 628
1.778 | cB211.072.64 | (64,128) | 137 | 374 | 2.1300(16)
cB211.14.64 | (64,128) | 190 | 437
cB211.25.48 | (48,96) | 253 | 314
cB211.25.32 | (32,64) | 253 | 400
1.836 | cC211.06.80 | (80,160) | 134 | 401 | 2.5045(17)
| cC211.20.48 | (48,96) | 246 | 890 |

Table 1: Simulation details for the $N_{f}=2+1+1$ ensembles by ETMC.

##  3 Determination of pseudoscalar meson decay constants with tmQCD

In the maximal tm (Mtm) formulation of LQCD the computation of pseudoscalar
decay constants does not require any (re)normalisation constant thanks to the
existence of a conserved current [1, 18]. It is thus sufficient to employ
correlation functions of the type
$C_{PP}(t)=(1/L^{3})\sum_{\vec{x},\vec{y}}\langle
0|P_{ff^{\prime}}(\vec{x},t)P^{\dagger}_{ff^{\prime}}(\vec{y},0)|0\rangle,$
where
$P_{ff^{\prime}}(x)=\bar{\chi}_{f}^{val}(x)\gamma_{5}\chi_{f^{\prime}}^{val}(x)$
(with flavours $\\{f,f^{\prime}\\}=\\{\ell,s,c\\}$) is the pseudoscalar
density operator. Then the leptonic decay constant of a PS-meson with mass
denoted by $M_{ps(ff^{\prime})}$ made out of valence quark flavours with bare
masses $\mu_{f}$ and $\mu_{f^{\prime}}$ is given by

$f_{ps}=(\mu_{f}+\mu_{f^{\prime}})\dfrac{\langle
0|P_{ff^{\prime}}|ps\rangle}{M_{ps(ff^{\prime})}\sinh(M_{ps(ff^{\prime})})}~{},$
(4)

which is automatically O$(a)$-improved owing to the maximal twist condition
($m_{0}^{val}=m_{0\ell}=m_{0h}=m_{crit}$).

In the present computation we make use of the quark mass results for the
$u/d$, $s$ and $c$ presented in Section V of Ref. [14] and corresponding to
the meson sector analysis. In this way we can employ the same mesonic
correlation functions as in Ref. [14] and hence determine correctly the error
propagation owing to the quark masses’ uncertainties. We recall that in the
framework of Mtm LQCD the renormalised quark mass of a quark flavour $f$ is
given by $m_{f}=\mu_{f}/Z_{P}$, where $Z_{P}$ is the renormalization constant
for the pseudoscalar density operator, the determination of which has been
presented in Refs [14, 19].

For the statistical and fit error analysis we have employed the jackknife
method. We have determined the continuum limit values at the physical light
quark mass $u/d$ by making use of simultaneous continuum and chiral fits. For
the estimation of the various sources of systematic uncertainty we repeat our
analysis by employing different kinds of fit ansätze regarding the chiral
extrapolation/interpolation to the physical light quark mass $u/d$.
Furthermore we perform several analyses by using data combinations
corresponding to two out of three lattice spacings and also by employing
different determinations for $Z_{P}$ that differ by O$(a^{2})$ effects. For a
given decay constant we thus obtain a distribution of results, with each
result corresponding to a different analysis and a total number of analyses
ranging from 32 to 96 (depending on the considered decay constant). From such
a distribution the mean value and the uncertainty of the final result are
estimated using the combination method and formulae discussed in Sec. V of
Ref. [14] (see there Eqs (38)–(43)).

##  4 Determinations of $f_{K}$ and $f_{K}/f_{\pi}$

In the determination of $f_{K}$ we first interpolate the decay constant
estimates to the strange quark mass and then we employ simultaneous continuum
and chiral fits of the quantity $f_{s\ell}$ against the light quark mass
$m_{\ell}$. We make use of two fit ansätze, namely the next-to-leading order
(NLO) SU(2) ChPT formula,
$f_{s\ell}=P_{0}\,(1-(3/4)\xi_{\ell}\log\xi_{\ell}+P_{1}\xi_{\ell}+P_{2}a^{2})\,\text{K}_{f_{K}}^{\text{FSE}}$
and a polynomial quadratic fit of the form
$f_{s\ell}=Q_{0}^{\prime}\,\left(1+Q_{1}^{\prime}m_{\ell}+Q_{2}^{\prime}m_{\ell}^{2}+Q_{3}^{\prime}a^{2}\right)\,\text{K}_{f_{K}}^{\text{FSE}}$.
In the first one it is set $\xi_{\ell}=(2B_{0}m_{\ell})/(4\pi f_{0})^{2}$
where $B_{0}$ and $f$ are the SU(2) ChPT low-energy constants (LECs) obtained
from the quark mass analysis, see Ref [14]. The factor
$\text{K}_{f_{K}}^{\text{FSE}}$ represents the estimation for the (small)
correction to our data due to finite size volume effects (FSE) following [20].

We work in a similar way for the determination of the ratio $f_{K}/f_{\pi}$
for which we make use of the following two fit ansätze:
$f_{s\ell}/f_{\ell\ell}=P_{0}^{\prime}\,(1+(5/4)\xi_{\ell}\log\xi_{\ell}+P_{1}^{\prime}\xi_{\ell}+P_{2}^{\prime}a^{2})\,\text{K}_{f_{K}/f_{\pi}}^{\text{FSE}}$
and
$f_{s\ell}/f_{\ell\ell}=Q_{0}^{\prime}\,\left(1+Q_{1}^{\prime}m_{\ell}+Q_{2}^{\prime}m_{\ell}^{2}+Q_{3}^{\prime}a^{2}\right)\,\text{K}_{f_{K}/f_{\pi}}^{\text{FSE}}$.
In Fig. 1 we show a representative plot out of several analyses concerning the
simultaneous chiral and continuum fits for $f_{K}$ (left panel) and
$f_{K}/f_{\pi}$ (right panel). In both cases the (NLO) SU(2) ChPT and
polynomial fits to our data are of good quality and the corresponding
continuum results at the physical point are compatible. After averaging over
results from all available analyses our preliminary results and the respective
error budgets for the two quantities are

$\displaystyle f_{K}^{\text{isoQCD}}$
$\displaystyle=155.3(0.9)_{\text{(stat+fit)}}(0.1)_{Z_{P}}(0.2)_{\text{chiral}}(1.4)_{\text{discr.}}(0.2)_{\text{FSE}}~{}[1.7]~{}\text{MeV}$
(5) $\displaystyle(f_{K}/f_{\pi})^{\text{isoQCD}}$
$\displaystyle=1.2023(38)_{\text{(stat+fit)}}(3)_{Z_{P}}(11)_{\text{chiral}}(8)_{\text{discr.}}(5)_{\text{FSE}}~{}[41],$
(6)

where separate errors in parentheses are due to the indicated sources of
uncertainty while the total error is shown in brackets. Notice that the
estimate for $(f_{K}/f_{\pi})^{\text{isoQCD}}$ has a total uncertainty of
about 0.34% and it shows nice agreement with the result of Ref. [9], namely
$(f_{K}/f_{\pi})^{\text{isoQCD}}=1.1995~{}(44)$, where an analysis of the same
data has been performed but in terms of the PS-meson masses. Note also that a
much more precise estimate for the $K$-decay constant, still in good agreement
with the result of Eq. (5), is obtained by

$f_{K}^{\text{isoQCD}}=(f_{K}/f_{\pi})^{\text{isoQCD}}\times
f_{\pi}^{\text{(isoQCD)}}=156.8(0.6)~{}\text{MeV}.$

Finally, by employing the estimate of the strong isospin effects correction
computed in Ref. [17] (in the GRS scheme [21]) we also obtain:

$f_{{K}^{\pm}}/f_{{\pi}^{\pm}}=1.1984~{}(41)~{}~{}\text{and}~{}~{}f_{{K}^{\pm}}=(f_{{K}^{\pm}}/f_{{\pi}^{\pm}})\times
f_{\pi}^{\text{(phys.)}}=156.3(0.6)~{}\text{MeV}$

|  |
---|---|---
Figure 1: Plots showing simultaneous chiral and continuum fits (continuum
limit curves are shown) for $f_{K}$ (left) and $f_{K}/f_{\pi}$ (right) against
the renormalised light quark mass
$m_{\ell}(\overline{\text{MS}},2~{}\text{GeV})$. Ansätze for chiral and
polynomial fits are given in the text. In both plots the vertical dotted line
indicates the physical value of $m_{u/d}$.

##  5 Determinations of $f_{D_{s}}$, $f_{D}$ and $f_{D_{s}}/f_{D}$

We compute $f_{D_{s}}$ by first interpolating the decay constant estimates to
the strange and charm quark masses before we employ a combined chiral and
continuum fit of the data for $f_{sc}$ in terms of $m_{\ell}$ and $a^{2}$. We
try two kinds of intermediate scaling variables that are the gradient flow
$w_{0}$ and the pseudocalar mass $M_{sc}$, the latter computed at each value
of $m_{\ell}$. The fit ansatz for both scaling variable choices is of the form
$f_{sc}=c_{0}\,\left(1+c_{1}m_{\ell}+c_{3}a^{2}\right)$ that is linear in
$m_{\ell}$ and it describes nicely our data as it can be appreciated by the
two plots in Fig. 2. Dependence on $m_{\ell}$, according to the expectations,
is quite weak. It should be also added that the continuum limit results from
both ways of analysis are in perfect agreement, nevertheless when $M_{sc}$ is
employed as the scaling variable cut off effects are clearly suppressed.

|  |
---|---|---
Figure 2: Plots showing simultaneous chiral and continuum fit for $f_{D_{s}}$
in terms of intermediate scaling variable $w_{0}$ (left) and of a PS-meson
mass $M_{sc}$ (right) against the renormalised light quark mass
$m_{\ell}(\overline{\text{MS}},2~{}\text{GeV})$. In both plots the vertical
dotted line indicates the physical value of $m_{u/d}$ and “C.L.” is for the
continuum limit curve.

For the calculation of the ratio $f_{D_{s}}/f_{D}$ we employ the following
three fit ansätze, where the first two ansätze are polynomial fits (linear or
quadratic) in $m_{\ell}$, namely $f_{sc}/f_{\ell
c}=\tilde{Q}_{0}(1+\tilde{Q}_{1}m_{\ell}+[\tilde{Q}_{2}m_{\ell}^{2}]+\tilde{Q}_{3}a^{2})$,
and the third one is based on the HMChPT prediction and takes the form
$f_{sc}/f_{\ell
c}=\tilde{P}_{0}(1+\frac{3}{4}(1+3\hat{g}^{2})\xi_{\ell}\log(\xi_{\ell})+\tilde{P}_{1}m_{\ell}+\tilde{P}_{2}a^{2})$,
where $\hat{g}=0.61(7)$ is obtained from the experimental measurement of the
$g_{D^{\star}D\pi}$. In Fig. 3 we present a representative analysis plot for
the ratio $f_{D{s}}/f_{D}$ against $m_{\ell}$ where all three kinds of
combined continnum and chiral fits are shown. We find that the HMCHPT fit
describes rather poorly our data. Therefore having data at (or close to) the
physical point we trust as for the central value and the error only fits of
good quality which in the present case are the two kinds of the polynomial
chiral fit ansatz.

---
Figure 3: Simultaneous continuum and chiral fits (continuum limit curves are
shown) for the ratio $f_{D_{s}}/f_{D}$ against the renormalised light quark
mass $m_{\ell}(\overline{\text{MS}},2~{}\text{GeV})$ using various chiral fit
ansätze. Vertical dotted line indicates the physical value of $m_{u/d}$.

Finally, we can compute $f_{D}$ directly and indirectly. In the first way, the
direct one, we follow an analysis similar to the $f_{D{s}}$ case i.e. we
employ two kinds of intermediate scaling variable. As for the chiral fit
ansatz we employ both polynomial and HMChPT fits. We observe that, similarly
to the $f_{D_{s}}$ case, the use of a PS-meson mass of the type $M_{\ell c}$
playing the role of intermediate scaling variable leads to suppressed
discretisation effects. Moreover as it happens for the analysis of the ratio
$f_{D{s}}/f_{D}$, also in the case of $f_{D}$ the HMChPT fit ansatz does not
provide a satisfactory fit quality. The indirect way for the computation of
$f_{D}$ consists in combining our results for the ratio and the $D_{s}$ decay
constants as follows $f_{D}=f_{D_{s}}/(f_{D_{s}}/f_{D})$.

Our preliminary results and error budget read:

$\displaystyle f_{D_{s}}$
$\displaystyle=248.9~{}(1.6)_{\text{(stat+fit)}}(0.5)_{Z_{P}}(0.2)_{\text{chiral}}(1.0)_{\text{discr.}}~{}[2.0]~{}\text{MeV}$
(7) $\displaystyle f_{D}$
$\displaystyle=210.1~{}(2.2)_{\text{(stat+fit)}}(0.1)_{Z_{P}}(0.4)_{\text{chiral}}(0.8)_{\text{discr.}}~{}[2.4]~{}\text{MeV}$
(8) $\displaystyle f_{D_{s}}/f_{D}$
$\displaystyle=1.1838~{}(90)_{\text{(stat+fit)}}(25)_{Z_{P}}(38)_{\text{chiral}}(57)_{\text{discr.}}~{}[115],$
(9)

where the total error for each quantity is shown in brackets. Notice that the
total relative errors for $f_{Ds}$, $f_{D}$ and $f_{D_{s}}/f_{D}$ are 0.8%,
1.1% and 1.0%, respectively.

##  6 Summary and results comparisons

In Table 2 we provide the comparison of the present work results (ETMC 21)
with older ETMC results (ETMC 14 [22]), the latter obtained with $N_{f}=2+1+1$
simulations but far from the physical point, and also with the FLAG 19
averages [11]. We would like to stress the higher precision of the ETMC 21
results with respect to the corresponding ETMC 14 ones owing to a remarkable
reduction of both statistical and systematic uncertainties. We also notice
that the ETMC 21 results compare well with the corresponding FLAG 19 averages.

Quantity | ETMC 21 | ETMC 14 | FLAG 19 | FLAG 19
---|---|---|---|---
| ($N_{f}=2+1+1$) | ($N_{f}=$ 2+1+1) | ($N_{f}=$ 2+1+1) | ($N_{f}=$ 2+1)
$(f_{K}/f_{\pi})^{\text{isoQCD}}$ | 1.2023(41) | 1.188(15) | - |
$f_{{K}^{\pm}}/f_{{\pi}^{\pm}}$ | 1.1984(41) | 1.184(16) | 1.1932(19) | 1.1917(37)
$f_{K}^{\text{isoQCD}}$ (MeV) | 155.3(1.7) | 155.0(1.9) | - | -
$f_{K}^{\text{isoQCD}}=(f_{K}/f_{\pi})^{\text{isoQCD}}\times f_{\pi}^{\text{(isoQCD)}}$ (MeV) | 156.8(0.6) | 154.9(1.9) | - | -
$f_{{K}^{\pm}}=(f_{{K}^{\pm}}/f_{{\pi}^{\pm}})\times f_{\pi}^{\text{(phys.)}}$ (MeV) | 156.3(0.6) | 154.4(2.0) | 155.7(0.3) | 155.7(0.7)
$f_{D_{s}}$ (MeV) | 248.9(2.0) | 247.2(4.1) | 249.9(0.5) | 248.0(1.6)
$f_{D_{s}}/f_{D}$ | 1.1838(115) | 1.192(22) | 1.1783(16) | 1.1740(70)
$f_{D}$ (MeV) | 210.1(2.4) | 207.4(3.8) | 212.0(0.7) | 209.0(2.4)
$(\dfrac{f_{D_{s}}}{f_{D}})/(\dfrac{f_{K}}{f_{\pi}})$ | 0.995(13) | 1.003(14) | - | -

Table 2: Comparison of (preliminary) results for the PS-meson decay constants
of the present work (ETMC 21) with previous ETMC results (ETMC 14 [22]) and
FLAG 19 averages [11].

Finally, by combining our results for the decay constants with the relevant
experimental inputs we provide estimates for several of the first and second
row elements of the CKM matrix. We get the following preliminary results:
$|V_{us}/V_{ud}|=0.2303(8)_{\text{th}}(3)_{\text{expt}}[8]$,
$|V_{us}|=0.2242(8)_{\text{th}}(3)_{\text{expt}}[8]$,
$|V_{cd}|=0.2199(25)_{\text{th}}(57)_{\text{expt}}[62]$ and
$|V_{cs}|=0.9871(79)_{\text{th}}(185)_{\text{expt}}[201]$. Thanks to the above
estimates for the unitarity checks of the first and the second CKM rows we
get:

$\displaystyle|V_{ud}|^{2}+|V_{us}|^{2}+|V_{ub}|^{2}-1$
$\displaystyle=-1.56(0.34)_{\text{th}}(0.62)_{\text{expt}}[0.71]\times
10^{-3}$ (10) $\displaystyle|V_{cd}|^{2}+|V_{cs}|^{2}+|V_{cb}|^{2}-1$
$\displaystyle=+2.3(1.6)_{\text{th}}(3.7)_{\text{expt}}[4.0]\times 10^{-2}.$
(11)

where $|V_{ub}|^{2}$ and $|V_{cb}|^{2}$ being of order $10^{-6}$ and
$10^{-4}$, respectively, have negligible impact to the present accuracy. Our
results lead to about 2$\sigma$ tension for the unitarity check of the first
row (at the per mille level) while they confirm the second row unitarity of
the CKM matrix at the percent level.

## Acknowledgments

We acknowledge PRACE (Partnership for Advanced Computing in Europe) for
awarding us access to the high-performance computing system Marconi and
Marconi100 at CINECA (Consorzio Interuniversitario per il Calcolo Automatico
dell’Italia Nord-orientale) under the grants Pra17-4394, Pra20-5171 and
Pra22-5171, and CINECA for providing us CPU time under the specific initiative
INFN-LQCD123. P.D. acknowledges support form the European Unions Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No. 813942 (EuroPLEx) and from INFN under the research project INFN-
QCDLAT. R.F. acknowledges the University of Rome Tor Vergata for the support
granted to the project PLNUGAMMA.

## References

  * [1] R. Frezzotti and G.C. Rossi, _Chirally improving wilson fermions. i: O(a) improvement_ , _JHEP_ 08 (2004) 007 [hep-lat/0306014].
  * [2] R. Frezzotti and G.C. Rossi, _Twisted-mass lattice QCD with mass non-degenerate quarks_ , _Nucl. Phys. Proc. Suppl._ 128 (2004) 193 [hep-lat/0311008].
  * [3] [ETMC] C. Alexandrou et al., _Simulating twisted mass fermions at physical light, strange and charm quark masses_ , _Phys. Rev. D_ 98 (2018) 054518 [1807.00495].
  * [4] [ETMC] A. Abdel-Rehim et al., _First physics results at the physical pion mass from $N_{f}=2$ Wilson twisted mass fermions at maximal twist_, _Phys. Rev. D_ 95 (2017) 094515 [1507.05068].
  * [5] [ALPHA] P. Dimopoulos, H. Simma and A. Vladikas, _Quenched B(K)-parameter from Osterwalder-Seiler tmQCD quarks and mass-splitting discretization effects_ , _JHEP_ 0907 (2009) 007 [0902.1074].
  * [6] D. Becirevic et al., _Exploring twisted mass lattice qcd with the clover term_ , _Phys. Rev._ D74 (2006) 034501 [hep-lat/0605006].
  * [7] Y. Iwasaki, _Renormalization Group Analysis of Lattice Theories and Improved Lattice Action: Two-Dimensional Nonlinear O(N) Sigma Model_ , _Nucl.Phys._ B258 (1985) 141.
  * [8] S. Aoki, R. Frezzotti and P. Weisz, _Computation of the improvement coefficient c(SW) to one loop with improved gluon actions_ , _Nucl. Phys. B_ 540 (1999) 501 [hep-lat/9808007].
  * [9] [ETMC] C. Alexandrou et al., _Ratio of kaon and pion leptonic decay constants with $N_{f}=2+1+1$ Wilson-clover twisted-mass fermions_, 2104.06747.
  * [10] [ETMC] J. Finkenrath et al., _Twisted mass gauge ensembles at physical values of the light, strange and charm quark masses_ , _PoS(LATTICE2021)284_ (2021) .
  * [11] [FLAG] S. Aoki et al., _FLAG Review 2019: Flavour Lattice Averaging Group (FLAG)_ , _Eur. Phys. J. C_ 80 (2020) 113 [1902.08191].
  * [12] K. Osterwalder and E. Seiler, _Gauge Field Theories on the Lattice_ , _Ann. Phys._ 110 (1978) 440.
  * [13] R. Frezzotti and G.C. Rossi, _Chirally improving Wilson fermions. II: Four-quark operators_ , _JHEP_ 10 (2004) 070 [hep-lat/0407002].
  * [14] [ETMC] C. Alexandrou et al., _Quark masses using twisted mass fermion gauge ensembles_ , 2104.13408.
  * [15] [ETMC] B. Kostrzewa et al., _Gradient-flow scale setting with $N_{f}=2+1+1$ Wilson-clover twisted-mass fermions_, _PoS(LATTICE2021)131_ (2021) .
  * [16] [FLAG] S. Aoki et al., _Review of lattice results concerning low-energy particle physics_ , _Eur. Phys. J. C_ 77 (2017) 112 [1607.00299].
  * [17] M. Di Carlo, D. Giusti, V. Lubicz, G. Martinelli, C.T. Sachrajda, F. Sanfilippo et al., _Light-meson leptonic decay rates in lattice QCD+QED_ , _Phys. Rev. D_ 100 (2019) 034514 [1904.08731].
  * [18] R. Frezzotti, P.A. Grassi, S. Sint and P. Weisz, _Lattice QCD with a chirally twisted mass term_ , _JHEP_ 0108 (2001) 058 [hep-lat/0101001].
  * [19] [ETMC] M. Di Carlo et al., _Renormalization constants of quark bilinear operators in QCD with dynamical up, down, strange and charm quarks_ , _PoS(LATTICE2021)399_.
  * [20] G. Colangelo, S. Durr and C. Haefeli, _Finite volume effects for meson masses and decay constants_ , _Nucl. Phys. B_ 721 (2005) 136 [hep-lat/0503014].
  * [21] J. Gasser, A. Rusetsky and I. Scimemi _Electromagnetic corrections in hadronic processes_ , _Eur. Phys. J._ C32 (2003) 97 [hep-ph/0305260 [hep-ph]].
  * [22] [ETMC] N. Carrasco et al., _Leptonic decay constants $f_{K},f_{D},$ and $f_{{D}_{s}}$ with $N_{f}=2+1+1$ twisted-mass lattice QCD_, _Phys. Rev._ D91 (2015) 054507 [1411.7908].