11institutetext: Observatoire de Genève, Université de Genève, Chemin Pegasi
51b, 1290 Sauverny, Switzerland 22institutetext: Departamento de Astronomía,
Universidad de Chile, Casilla 36-D, Santiago, Chile 33institutetext:
Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, CV4
7AL, UK 44institutetext: Centre for Exoplanets and Habitability, University of
Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK 55institutetext: Institute of
Planetary Research, German Aerospace Center, Rutherfordstrasse 2, 12489
Berlin, Germany 66institutetext: School of Physics and Astronomy, University
of Leicester, LE1 7RH, UK 77institutetext: Center for Astronomy and
Astrophysics, TU Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany
88institutetext: Astronomy Unit, Queen Mary University of London, Mile End
Road, London E1 4NS, UK 99institutetext: Astrophysics Group, Cavendish
Laboratory, J.J. Thomson Avenue, Cambridge CB3 0HE, UK 1010institutetext:
Department of Physics, Kavli Institute for Astrophysics and Space Research,
MIT, 77 Mass. Ave, Cambridge, MA 02139, USA 1111institutetext: Centro de
Astrofísica y Tecnologías Afines (CATA), Casilla 36-D, Santiago, Chile
1212institutetext: Facultad de Ingeniería y Ciencias, Universidad Adolfo
Ibáñez, Av.Diagonal las Torres 2640, Peñalolén, Santiago, Chile
1313institutetext: Millennium Institute for Astrophysics, Chile
1414institutetext: Instituto de Astronomía, Universidad Católica del Norte,
Angamos 0610, 1270709, Antofagasta, Chile

# NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star

Nolan Grieves , NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star
<EMAIL_ADDRESS>Louise D. Nielsen NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant star Jose I. Vines NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant star Edward M. Bryant NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Samuel Gill NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star
François Bouchy NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Monika
Lendl NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Daniel
Bayliss NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Philipp
Eigmueller NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Damien
Segransan NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Jack S.
Acton NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star David
R. Anderson NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Matthew
R. Burleigh NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Sarah
L. Casewell NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star
Alexander Chaushev NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Benjamin F. Cooke NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Edward
Gillen Winton FellowNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Michael
R. Goad NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star
Maximilian N. Günther Juan Carlos Torres FellowNGTS-13b: A hot 4.8 Jupiter-
mass planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant star Beth A. Henderson NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant star Aleisha Hogan NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant star James S. Jenkins NGTS-13b: A hot 4.8 Jupiter-mass
planet transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet
transiting a subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Douglas R. Alves NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Andrés Jordán NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star James
McCormac NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star
Maximiliano Moyano NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Didier Queloz NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Liam Raynard NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Julia V. Seidel NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Alexis M. S. Smith NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Rosanna H. Tilbrook NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Stephane Udry NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
star Richard G. West NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a
subgiant starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star Peter
J. Wheatley NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant
starNGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star

(Received 2 October 2020; accepted 6 January 2020)

We report the discovery of the massive hot Jupiter NGTS-13b by the Next
Generation Transit Survey (NGTS). The V = 12.7 host star is likely in the
subgiant evolutionary phase with $\log g_{*}$ = 4.04 $\pm$ 0.05, $T_{\rm eff}$
= 5819 $\pm$ 73 K, $M_{*}$ = 1.30${}^{+0.11}_{-0.18}$ $\mathrm{M_{\odot}}$,
and $R_{*}$ = 1.79 $\pm$ 0.06 $\mathrm{R_{\odot}}$. NGTS detected a transiting
planet with a period of $P$ = 4.12 days around the star, which was later
validated with the Transiting Exoplanet Survey Satellite (TESS; TIC
454069765). We confirm the planet using radial velocities from the CORALIE
spectrograph. Using NGTS and TESS full-frame image photometry combined with
CORALIE radial velocities we determine NGTS-13b to have a radius of $R_{P}$ =
1.142 $\pm$ 0.046 $\mathrm{R_{\rm Jup}}$, mass of $M_{P}$ = 4.84 $\pm$ 0.44
$\mathrm{M_{\rm Jup}}$, and eccentricity $e$ = 0.086 $\pm$ 0.034. Some
previous studies suggest that $\sim$4 $\mathrm{M_{\rm Jup}}$ may be a border
between two separate formation scenarios (e.g., core accretion and disk
instability) and that massive giant planets share similar formation mechanisms
as lower-mass brown dwarfs. NGTS-13b is just above 4 $\mathrm{M_{\rm Jup}}$
making it an important addition to the statistical sample needed to understand
the differences between various classes of substellar companions. The high
metallicity, $\left[{\rm Fe}/{\rm H}\right]$ = 0.25 $\pm$ 0.17, of NGTS-13
does not support previous suggestions that massive giants are found
preferentially around lower metallicity host stars, but NGTS-13b does support
findings that more massive and evolved hosts may have a higher occurrence of
close-in massive planets than lower-mass unevolved stars.

###### Key Words.:

Planets and satellites: detection – Planets and satellites: individual:
(NGTS-13b) – techniques: photometric – techniques: radial velocities

##  1 Introduction

The first extrasolar planets discovered were a new class of planets not found
in our solar system: hot Jupiters – Jupiter-mass planets with orbital periods
$\leq$ 10 days, e.g., Mayor & Queloz (1995). Planet formation and evolution
theories were abruptly altered to account for these distinct systems, e.g.,
Boss (1996). However, there is still no consensus on the predominant formation
mechanism of these massive close-in planets. Currently three main hypotheses
are suggested as possible mechanisms for how hot Jupiters could have formed
and now reside so close to their host star: in situ formation, disk migration,
and high-eccentricity tidal migration. For a review of these scenarios see
Dawson & Johnson (2018).

There are two mechanisms proposed for giant planet formation in general:
gravitational instability – through its own self-gravity the protoplanetary
disk fragments into giant gaseous protoplanets that contract and collapse to
form giant planets (e.g., Boss, 1997), and core accretion – planetesimal
collisions form rocky protoplanet cores that accrete gas from the
protoplanetary disk (e.g., Pollack et al., 1996). In situ gravitational
instability is not plausible at short orbital periods because the disk’s fast
rotation and high temperature (thermal pressure) support the local gas against
gravitational collapse (Dawson & Johnson, 2018). Core accretion can only occur
if there is enough mass in the feeding zone to grow a sufficiently large core
($\sim$10 $\mathrm{M_{\oplus}}$; Rafikov, 2006; Lee et al., 2014; Piso et al.,
2015) within the disk’s lifetime, but at short orbital periods the feeding
zone is smaller (smaller Hill radii) making it implausible to form a massive
core in normal disk conditions (Lee & Chiang, 2016). However, in situ core
accretion of hot Jupiters may be possible with highly efficient radial drift
of solids to the inner disk (Dawson & Johnson, 2018).

A more common line of thought is that hot Jupiters formed at larger orbital
periods and migrated inward. One mechanism for this migration is via disk
migration by exchanging angular momentum with the disk through corotation and
Lindblad torques (e.g., Lin et al., 1996). Another proposed mechanism for
migration is high-eccentricity tidal migration, where a giant planet has
orbital angular momentum extracted from it by a perturber, pushing the planet
to a highly elliptical orbit, then the planet reduces its orbital energy
through tidal interactions with the central star causing it to circularize
into a close-in orbit (e.g., Rasio & Ford, 1996).

These various formation theories should have differing effects on the
characteristics of hot Jupiters. Slow and smooth migration in the
circumstellar disk from simple angular momentum loss should produce aligned
and circularized orbits (D’Angelo et al., 2003), but more dynamic migration,
such as high-eccentricity tidal migration, should produce planets with higher
eccentricites in more misaligned orbits (e.g., Triaud et al., 2010). However,
for hot Jupiters with orbital periods shorter than 5 days the eccentricity
will be erased in a few billion years due to tidal circularization (Adams &
Laughlin, 2006). Orbital misalignment, a difference in a star’s rotation axis
and a planet’s orbital axis (obliquity), evolves more slowly than eccentricity
(Hut, 1981; Winn et al., 2005; Barker & Ogilvie, 2009) and therefore can be
detected in relatively older stars. There are various methods for measuring
obliquity, but particularly for transiting planets radial velocity (RV)
observations can be taken during the transit to find the spin-orbit alignment
via the Rossiter-McLaughin effect (Rossiter, 1924; McLaughlin, 1924).

The parameter distribution of the entire sample of hot Jupiters such as their
semi-major axes and ages also give insight into which formation mechanisms
dominate in general and in specific environments. It is thus important to
continue to build a sample of planets (and host stars) with well-characterized
parameters. When a planet is detected by both the transit and RV methods
multiple characteristics of the planet can be well-constrained compared to
only detecting it through one method. Thus, large photometric surveys are
valuable as they can view large samples of stars and identify those with
transiting planets, which can then be followed up with RV instruments that
typically can only observe one star at a time.

The Next Generation Transit Survey (NGTS; Wheatley et al., 2018) is a ground-
based photometric survey located in Chile at the European Southern Observatory
(ESO) Paranal site. NGTS consists of twelve fully robotic telescopes each with
a 20 cm photometric aperture and a wide 8 deg2 field-of-view making the survey
capable of observing a large sample of stars. Here we present the discovery of
a massive hot Jupiter orbiting a subgiant star. This planet was first detected
by NGTS, subsequently verified by the space-based all-sky Transiting Exoplanet
Survey Satellite (TESS; Ricker et al., 2015), and confirmed through follow-up
RV measurements with the CORALIE (Queloz et al., 2001) spectrograph. In
Section 2 we present our observations, in Section 3 we describe our analysis
and present our results, in Section 4 we discuss our results, and we present
our conclusions in Section 5.

##  2 Observations

Figure 1: TESS Full-Frame Image cutout (15 $\times$ 15 pixels) in the region
of NGTS-13 (TIC 454069765; T=12.11) from which the TESS light curve was
generated. The red boxes indicate the target pixel aperture mask used to
extract the photometry and the white boxes without a border are the pixels
used to determine the sky background. The red star indicates the position of
NGTS-13. The yellow crosses indicate TIC 454069767 ($\Delta$T=1.68, 22.5”
west) and TIC 454069772 ($\Delta$T=0.6, 57” southeast). Blue pluses denote the
positions of stars in $GAIA$ DR2 within 1’ of NGTS-13.

Figure 2: Top: NGTS observations (binned to 30-minutes for visual purposes)
with the EXOFASTv2 model (see Section 3.2) in red. Middle: TESS Full Frame
Image (FFI) observations after the detrending described in Section 2.2 was
applied with the the EXOFASTv2 model in red. Bottom: NGTS binned to 30-minutes
and TESS FFI photometry phased to the 4.12 day period of NGTS-13b. The NGTS
data are binned to 30-minutes and then phased to the period. The red lines
show the EXOFASTv2 model.

###  2.1 NGTS Photometry

NGTS-13 was observed using a single NGTS telescope from UT 2017 December 10 to
UT 2018 August 7. During this period 255,709 observations were taken with
exposure times of 10 seconds. Image reduction and photometry were performed
with the NGTS pipeline (Wheatley et al., 2018) using CASUtools
software111http://casu.ast.cam.ac.uk/surveys-projects/software-release. The
lightcurve was detrended and systematic effects were removed using the SysRem
algorithm (Tamuz et al., 2005).

The NGTS lightcurves were searched for periodic transit-like signals using an
implementation of the box-fitting least-squares (BLS) algorithm (Kovács et
al., 2002). We detected a strong peak in the BLS periodogram at 4.12 days and
folding the NGTS photometry on this period revealed a clear transit signal.

We binned the NGTS data to 2-minute bins for our analysis as described in
Section 3.2. NGTS-13 has some nearby companions (see Figure 1) and the NGTS
photometry is expected to suffer from some dilution, which we discuss in
Section 3.2.2. The NGTS photometry binned to 30 minutes for display purposes
is shown in Figure 2. A full table of reduced NGTS photometric observations is
available in a machine-readable format from the online journal.

###  2.2 TESS Photometry

TESS observed NGTS-13 (TIC 454069765) in Sector 10 from UT 2019 March 28 to UT
2019 April 21 with 30-minute cadence Full Frame Images (FFIs). We created a
lightcurve via aperture photometry of the TESS 30-minute cadence FFIs
following the method set out in Gill et al. (2020). We determined a threshold
for target pixels and for background pixels based on an iterative sigma-
clipping to determine the median and standard deviation of the background
pixels. To exclude the neighboring T = 12.8 star (TIC 454069772, 57”
southeast, see Figure 1), we only include pixels in our final aperture for
which neighbouring pixels closer to the center of our target show a higher
illumination. For the final pixel mask of NGTS-13, 7 pixels were selected (see
Figure 1). To remove the systematic trends from the TESS light curve, we mask
out the transits, and then fit a linear spline interpolating across the
positions of the transits. Before removing systematic trends, we estimate a
signal-to-noise ratio (SNR) as the transit depth relative to the standard
deviation of out-of-transit observations and find a single TESS transit
SNR$\sim$5.4, which corresponds to SNR$\sim$13.3 when considering all 6
transits in the TESS data (SNR $\propto$ $\sqrt{N_{\mathrm{T}}}$). The TESS
photometry is displayed in Figure 2. A full table of reduced TESS photometric
observations is available in a machine-readable format from the online
journal.

###  2.3 CORALIE Spectroscopy

We obtained 13 observations of NGTS-13 with the high resolution echelle
spectrograph CORALIE on the Swiss 1.2 m Euler telescope (Queloz et al., 2001)
at La Silla Observatory, Chile between UT 2019 December 12 and UT 2020 January
23. CORALIE has a resolution of $R$ $\sim$ 60,000 and is fed by two fibers: a
2 arcsec on-sky science fiber encompassing the star and another fiber that can
either connect to a Fabry-Pérot etalon for simultaneous wavelength calibration
or on-sky for background subtraction of sky flux. We observed NGTS-13 in the
simultaneous Fabry-Pérot wavelength calibration mode using an exposure time of
2700 seconds. The spectra were reduced with the CORALIE standard reduction
pipeline and RVs were computed for each epoch by cross-correlating with a
binary G2 mask (Pepe et al., 2002). We co-added the 13 CORALIE spectroscopic
observations by correcting each epoch to the stellar rest frame with its
corresponding barycentric RV to create a single high signal-to-noise spectrum.
We list the CORALIE RVs in Table 1. The RVs phased to the 4.12 day period are
displayed in Figure 3.

##  3 Analysis

Figure 3: CORALIE Radial Velocities for NGTS-13 phased to the 4.12 day period of NGTS-13b subtracted by the relative RV offset $\gamma_{\rm rel}$ (Table 4). The red line shows the Keplerian fit for the planet from EXOFASTv2 (see Section 3.2). BJD TDB | RV [m s-1] | RV error [m s-1]
---|---|---
2458829.841574 | 26749.44 | 51.75
2458839.806533 | 27431.37 | 59.75
2458841.801379 | 26588.94 | 71.16
2458855.811893 | 27554.10 | 85.54
2458857.832753 | 26497.28 | 54.13
2458858.814847 | 26890.46 | 65.16
2458860.806340 | 27190.32 | 59.13
2458862.815160 | 26824.62 | 51.75
2458863.794625 | 27457.32 | 64.68
2458864.726425 | 27354.61 | 78.83
2458868.781343 | 27378.28 | 68.78
2458869.733432 | 26669.54 | 88.67
2458871.821928 | 27459.18 | 60.68
Table 1: CORALIE Radial Velocities for NGTS-13. Time is in Barycentric Julian Date / Barycentric Dynamical Time. $T_{\rm eff}$ | $\left[{\rm Fe}/{\rm H}\right]$ | $\log g_{*}$ | Source
---|---|---|---
K | dex | dex |
5800 $\pm$ 80 | 0.20 $\pm$ 0.08 | 4.1 $\pm$ 0.2 | Wavelet
5714 $\pm$ 110 | 0.43 $\pm$ 0.09 | 4.11 $\pm$ 0.12 | SpecMatch
5757 $\pm$ 153 | 0.315 $\pm$ 0.205 | 4.105 $\pm$ 0.205 | Combined
Table 2: Table of stellar atmospheric parameters derived from CORALIE spectra and combining the analysis results of both SpecMatch-Emp (Yee et al., 2017) and the wavelet analysis described in Gill et al. (2018, 2019). We use the combined results as input priors into our global analysis. Parameter | Description (Units) | Value | Source
---|---|---|---
Astrometric Properties |  |  |
$\alpha$ | Right Ascension | 11:44:57.66 |  $Gaia$ DR2
$\delta$ | Declination | -38:08:22.98 |  $Gaia$ DR2
$\mu_{\text{R.A.}}$ | R.A. proper motion (mas y-1) | -21.126 $\pm$ 0.046 |  $Gaia$ DR2
$\mu_{\text{DEC}}$ | DEC proper motion (mas y-1) | -0.988 $\pm$ 0.037 |  $Gaia$ DR2
$\varpi$ | Parallax (mas) | 1.516 $\pm$ 0.035 |  $Gaia$ DR2
Photometric Properties |  |  |
V | (mag) | 12.704 $\pm$ 0.021 | APASS
B | (mag) | 13.419 $\pm$ 0.010 | APASS
J | (mag) | 11.481 $\pm$ 0.020 | 2MASS
H | (mag) | 11.244 $\pm$ 0.030 | 2MASS
K | (mag) | 11.141 $\pm$ 0.020 | 2MASS
SDSS g | (mag) | 13.006 $\pm$ 0.022 | SDSS
SDSS i | (mag) | 12.363 $\pm$ 0.027 | SDSS
SDSS r | (mag) | 12.514 $\pm$ 0.065 | SDSS
W1 | (mag) | 11.095 $\pm$ 0.030 | WISE
W2 | (mag) | 11.148 $\pm$ 0.030 | WISE
W3 | (mag) | 11.125 $\pm$ 0.107 | WISE
NUV | (mag) | 18.471 $\pm$ 0.066 | GALEX
Table 3: NGTS-13 (2MASS J11445767-38082292 / TIC 454069765) astrometric and
photometric properties. The astrometric and photometric properties as well as
the spectroscopically derived atmospheric parameters in Table 2 are used as
inputs into our EXOFASTv2 global fit including the SED fit (see section 3.2).

###  3.1 Spectral Analysis

Figure 4: The Spectral Energy distribution model (black line) of NGTS-13 from
the EXOFASTv2 fit with broad band averages (blue circles) and broad band
measurements (red). The error bars in wavelength denote the bandwidth of the
corresponding filter and the error bars in flux denote the measurement
uncertainty.

We used wavelet analysis to extract stellar atmospheric parameters from the
co-added 13 CORALIE spectroscopic observations of NGTS-13 using the
methodology set out in Gill et al. (2018, 2019). With the wavelet analysis we
find NGTS-13 to have surface rotational velocity $v\sin i_{*}$ = 6.2 $\pm$ 1.2
km s-1. We also derived stellar atmospheric parameters using SpecMatch-Emp
(Yee et al., 2017). SpecMatch-Emp uses a large library of stars with well-
determined parameters to match the input spectra and derive spectral
parameters. We use a spectral region that includes the Mg I b triplet (5100 -
5400 $\AA$) to match our spectra. SpecMatch-Emp uses $\chi^{2}$ minimisation
and a weighted linear combination of the five best matching spectra in the
SpecMatch-Emp library to determine $T_{\rm eff}$, $\log g_{*}$, and
$\left[{\rm Fe}/{\rm H}\right]$. We combine the wavelet analysis results and
SpecMatch-Emp results to create wide uncertainties that include the
uncertainty range of both methods to use as priors for their final derivation
as described below. We present these spectroscopically derived parameters from
both methods and the adopted combined values in table 2.

###  3.2 Global Analysis

We derive both planet and stellar parameters using EXOFASTv2 (Eastman et al.,
2013; Eastman, 2017; Eastman et al., 2019), which can globally fit all
photometry and RV data. A full description of EXOFASTv2 is given in Eastman et
al. (2019), which can fit any number of transits and RV sources while
exploring the vast parameter space through a differential evolution Markov
Chain coupled with a Metropolis-Hastings Monte Carlo sampler. A built-in
Gelman-Rubin statistic (Gelman & Rubin, 1992; Gelman et al., 2003; Ford, 2006)
is used to check the convergence of the chains.

We use Modules for Experiments in Stellar Astrophysics (MESA) Isochrones &
Stellar Tracks (MIST; Choi et al., 2016; Dotter, 2016) isochrones and the
spectral energy distribution (SED) within the EXOFASTv2 fit to determine host
star parameters. For the SED fit within EXOFASTv2 (Stassun & Torres, 2016) we
use photometry from APASS DR9 BV (Henden et al., 2016), 2MASS JHK (Skrutskie
et al., 2006), SDSS DR12 gri (Alam et al., 2015), ALL-WISE W1, W2 and W3
(Wright et al., 2010), and GALEX NUV (Bianchi et al., 2011), which are
presented in Table 3. The SED fit within EXOFASTv2 puts sytematic floors on
the broadband photometry errors (Stassun & Torres, 2016). The EXOFASTv2 SED
fit is presented in Figure 4. We use our spectroscopically derived $T_{\rm
eff}$, $\log g_{*}$, and $\left[{\rm Fe}/{\rm H}\right]$ values as priors in
the global model as well as a Gaussian prior for parallax from the value and
uncertainty in $Gaia$ DR2. We put an upper limit on the extinction based on
the maximum in line-of-sight value from the Schlafly & Finkbeiner (2011)
galaxy dust map. All other fitted and derived parameters from our EXOFASTv2
model have conservative physical boundaries that are detailed in Table 3 of
Eastman et al. (2019), which also gives a thorough explanation of each
parameter.

Due to the long exposure time of the TESS full frame exposures we must account
for smearing of the light curve. Therefore, within EXOFASTv2 we specify an
exposure time of 30 minutes and average over 10 data points to integrate a
model over the exposure time equivalent to a midpoint Riemann sum (Eastman et
al., 2019). We assign the TESS FFI photometry to the built-in EXOFASTv2 ”TESS”
band and the NGTS photometry to the ”R” band for the computation of the limb
darkening coefficients. We found the EXOFASTv2 ”R” band to have the most
similar bandpass to NGTS in the EXOFASTv2 framework, which only allows filters
defined in Claret & Bloemen (2011). We tested the use of this filter by
running a fit without TESS data and NGTS data assigned to the ”TESS” filter;
we obtained limb darkening parameters within uncertainties to those presented
in Table 4. We also ran a fit including TESS data assigned to the ”TESS”
filter and the NGTS data assigned to the ”Kepler” filter; we again found limb
darkening parameters within uncertainties to those presented. We account for
dilution of the photometry lightcurves, which we describe in Section 3.2.2.
For our EXOFASVTv2 global model and Markov chain Monte Carlo (MCMC) fit we
used 56 walkers (2 $\times$ n${}_{\text{parameters}}$), or chains, and allowed
the fit to run until convergence (132776 steps).

Parameter | Unit | Value | Value (High $\mathrm{M_{\odot}}$) | Value (Low $\mathrm{M_{\odot}}$)
---|---|---|---|---
|  | (Adopted Solution) | (60% Probability) | (40% Probability)
Stellar Parameters |  |  |  |
$M_{*}$ | Mass ($\mathrm{M_{\odot}}$) | $1.30^{+0.11}_{-0.18}$ | $1.34^{+0.05}_{-0.05}$ | $1.16^{+0.05}_{-0.07}$
$R_{*}$ | Radius ($\mathrm{R_{\odot}}$) | $1.788^{+0.057}_{-0.057}$ | $1.792^{+0.051}_{-0.049}$ | $1.787^{+0.052}_{-0.052}$
$L_{*}$ | Luminosity ($\mathrm{L_{\odot}}$) | $3.29^{+0.21}_{-0.21}$ | $3.34^{+0.19}_{-0.19}$ | $3.25^{+0.19}_{-0.20}$
$\rho_{*}$ | Density (cgs) | $0.31^{+0.04}_{-0.04}$ | $0.33^{+0.03}_{-0.03}$ | $0.28^{+0.03}_{-0.03}$
$\log g_{*}$ | Surface gravity (cgs) | $4.04^{+0.05}_{-0.05}$ | $4.06^{+0.02}_{-0.02}$ | $4.00^{+0.03}_{-0.03}$
$T_{\rm eff}$ | Effective Temperature (K) | $5819^{+73}_{-73}$ | $5831^{+62}_{-72}$ | $5799^{+64}_{-73}$
$\left[{\rm Fe}/{\rm H}\right]$ | Metallicity (dex) | $0.25^{+0.17}_{-0.17}$ | $0.28^{+0.13}_{-0.14}$ | $0.12^{+0.16}_{-0.18}$
Age | Age (Gyr) | $4.23^{+2.65}_{-1.59}$ | $3.68^{+0.58}_{-0.62}$ | $6.55^{+1.30}_{-0.89}$
EEP | Equal Evolutionary Point | $413^{+47}_{-18}$ | $404^{+9}_{-14}$ | $453^{+4}_{-5}$
AV | V-band extinction (mag) | $0.164^{+0.044}_{-0.044}$ | $0.159^{+0.036}_{-0.064}$ | $0.149^{+0.043}_{-0.071}$
$d$ | Distance (pc) | $657^{+15}_{-15}$ | $659^{+15}_{-14}$ | $655^{+15}_{-14}$
Planet Parameters |  |  |  |
$M_{P}$ | Mass ($\,M_{\rm J}$) | $4.84^{+0.44}_{-0.44}$ | $5.05^{+0.26}_{-0.26}$ | $4.54^{+0.26}_{-0.29}$
$R_{P}$ | Radius ($\,R_{\rm J}$) | $1.142^{+0.046}_{-0.046}$ | $1.142^{+0.036}_{-0.034}$ | $1.143^{+0.040}_{-0.038}$
$P$ | Period (days) | $4.119027^{+0.000023}_{-0.000023}$ | $4.119027^{+0.000020}_{-0.000020}$ | $4.119027^{+0.000020}_{-0.000020}$
$T_{C}$ | Time of conjunction ($\rm{BJD_{TDB}}$) | $2458572.8247^{+0.0018}_{-0.0018}$ | $2458572.8247^{+0.0015}_{-0.0015}$ | $2458572.8246^{+0.0015}_{-0.0015}$
$a$ | Semi-major axis (AU) | $0.0549^{+0.0015}_{-0.0025}$ | $0.0556^{+0.0007}_{-0.0007}$ | $0.0529^{+0.0007}_{-0.0011}$
$i$ | Inclination (Degrees) | $88.7^{+1.2}_{-1.2}$ | $88.8^{+0.8}_{-1.1}$ | $87.9^{+1.4}_{-1.3}$
$e$ | Eccentricity | $0.086^{+0.034}_{-0.034}$ | $0.091^{+0.023}_{-0.022}$ | $0.080^{+0.023}_{-0.022}$
$\omega_{*}$ | Argument of Periastron (Degrees) | $-146.4^{+24.2}_{-24.2}$ | $-143.0^{+13.0}_{-18.0}$ | $-161.0^{+24.0}_{-29.0}$
$K$ | RV semi-amplitude (m/s) | $526.1^{+29.5}_{-29.5}$ | $526.0^{+24.0}_{-24.0}$ | $525.0^{+25.0}_{-27.0}$
$\delta$ | Transit depth (fraction) | $0.0043^{+0.0004}_{-0.0004}$ | $0.0043^{+0.0001}_{-0.0001}$ | $0.0043^{+0.0001}_{-0.0001}$
$\tau$ | Ingress/egress transit duration (days) | $0.0150^{+0.0022}_{-0.0022}$ | $0.0137^{+0.0009}_{-0.0003}$ | $0.0143^{+0.0018}_{-0.0008}$
$T_{14}$ | Total transit duration (days) | $0.2176^{+0.0030}_{-0.0030}$ | $0.2173^{+0.0021}_{-0.0020}$ | $0.2181^{+0.0024}_{-0.0022}$
$b$ | Transit Impact parameter | $0.156^{+0.145}_{-0.145}$ | $0.150^{+0.130}_{-0.100}$ | $0.240^{+0.150}_{-0.160}$
$\rho_{P}$ | Density (cgs) | $4.02^{+0.55}_{-0.55}$ | $4.20^{+0.43}_{-0.41}$ | $3.77^{+0.46}_{-0.43}$
log gP | Surface gravity (cgs) | $3.96^{+0.06}_{-0.06}$ | $3.98^{+0.03}_{-0.03}$ | $3.94^{+0.04}_{-0.04}$
$T_{eq}$ | Equilibrium temperature (K) | $1605^{+30}_{-30}$ | $1595^{+21}_{-22}$ | $1627^{+24}_{-26}$
$S$ | Insolation Flux ($S_{\oplus}$) | $1096^{+90}_{-90}$ | $1064^{+58}_{-59}$ | $1153^{+70}_{-70}$

Wavelength Parameters |  | R (NGTS) | TESS
---|---|---|---
$u_{1}$ | linear limb-darkening coeff | $0.314^{+0.047}_{-0.047}$ | $0.331^{+0.049}_{-0.050}$
$u_{2}$ | quadratic limb-darkening coeff | $0.248^{+0.049}_{-0.049}$ | $0.307^{+0.049}_{-0.049}$

Transit Parameters |  | NGTS | TESS
---|---|---|---
$\sigma^{2}$ | Added Variance | $0.00000501^{+0.00000021}_{-0.00000020}$ | $0.00000034^{+0.00000003}_{-0.00000003}$
$F_{0}$ | Baseline flux | $1.000270^{+0.000047}_{-0.000047}$ | $0.999972^{+0.000026}_{-0.000026}$

RV Parameters |  | CORALIE
---|---|---
$\gamma_{\rm rel}$ | Relative RV Offset (m/s) | $27029^{+17}_{-17}$
$\sigma_{J}$ | RV Jitter (m/s) | $0^{+40}_{-0}$
$\sigma_{J}^{2}$ | RV Jitter Variance | $-800^{+2400}_{-1100}$
o$-$c rms | Root Mean Square of residuals (m s-1) | 40.6

Table 4: NGTS-13 stellar and companion parameters derived from our EXOFASTv2
modeling and analysis described in Section 3. We present two solution split
between the peaks of our bimodal stellar mass posterior distribution, as
described in Section 3.2.1. The final adopted solution is in bold and uses the
median value and upper and lower boundaries of the 68% confidence interval we
computed, as described in Section 3.2.1. For the wavelength, transit, and RV
parameters listed in the lower part of the table we found the two solutions to
be very similar and well below uncertainties and only present the median
values of all the MCMC steps. The boundaries (besides the specific priors
mentioned in Section 3) and thorough explanations of all fitted and derived
parameters from EXOFASTv2 are detailed in Table 3 of Eastman et al. (2019).

####  3.2.1 Bimodal Stellar Mass and Age

Figure 5: Stellar mass and age posterior distributions from EXOFASTv2. These
distributions show the relative probability between the two peaks of the
bimodal distributions, which is approximately 60% to 40% in favor of a more
massive, younger system. We split the stellar mass distribution in the valley
at 1.25 M⊙ to create two solutions presented in table 4. For our final
solution we adopt the median value (vertical black lines) and upper and lower
boundaries (dashed vertical lines) of the 68% confidence intervals we
computed, as described in Section 3.2.1.

We inspected the posterior distributions of each fitted and derived parameter
and found bimodal distributions in both the stellar mass and age. This type of
bimodality between two different stellar ages and masses occurs due to a
degeneracy between two MIST isochrone tracks that EXOFASTv2 interpolates
between and has been seen in several recent studies (e.g., Ikwut-Ukwa et al.,
2020; Pepper et al., 2020; Carmichael et al., 2020a, b). Similar to those
studies we split the host star mass posterior distribution at the valley
between the two peaks, or 1.25 M⊙, and extract two separate solutions, both of
which are presented in Table 3. We show the stellar mass and age posterior
distributions in Figure 5. The higher mass solution is moderately more
probable with a probability of 60% compared to 40% for the lower mass
solution.

However, for our final adopted solution instead of picking one peak we define
wider uncertainties by computing confidence intervals for each posterior
distribution based on the smallest area in the posterior distribution function
that contains 68% of the MCMC steps. We compute these confidence intervals
using a subroutine within the open source python package
MC3222https://mc3.readthedocs.io/en/latest/index.html (Cubillos et al., 2017).
In summary the code orders the posterior distribution by value and creates a
smoothed posterior density distribution (PDF), computes a cumulative
distribution function (CDF) of the ordered PDF, and computes where the CDF
reaches 68% of its maximum value. The minimum and maximum parameter values
inside this interval are the boundaries of the 1$\sigma$ confidence interval.
We use the median value within the confidence interval as our final value for
each parameter and the minimum and maximum values within the confidence
interval as their uncertainties, which are in bold in Table 4. Figure 5
displays these values for the host star mass and age.

####  3.2.2 Dilution

NGTS-13 has several close companions including two stars only slightly dimmer:
TIC 454069767 (TESS mag = 14.1, 22.5” west) and TIC 454069772 (TESS mag =
13.4, 57” south-east), see Figure 1. We expect the TESS photometry of NGTS-13
to be completely diluted by the companion 22.5” away and partially diluted by
the companion 57.3” away. One other close companion has a TESS mag $<$ 18,
while all other companions are dimmer and not expected to dilute at all.

To account for dilution we calculate a dilution factor, given as a contrast
ratio $C$, between the combined flux of all contaminating sources and the flux
of the target star:

$C=F_{cont}/F_{target}$ (1)

where $F_{target}$ and $F_{cont}$ are the flux contributions within the
photometric aperture from the target and contaminating stars respectively. We
follow the method detailed in Bryant et al. (2020) to estimate dilution.

We determine the point-spread function (PSF) in the region of the camera
surrounding NGTS-13 and extract 2D Gaussian shape parameters (Pál, 2009) using
the Cluster Difference Imaging Photometric Survey (CDIPS; Bouma et al., 2019).
For NGTS-13 we find shape parameters of $S$ = 1.88 $\pm$ 0.047 (FHWM = 1.71
pixels), $D$ = 0.033 $\pm$ 0.011, and $K$ = 0.169 $\pm$ 0.007. We consider all
the stars in TICv8 (Stassun et al., 2019) within 1’ of NGTS-13 and find a
dilution factor of 19.10${}^{+0.05}_{-0.04}$% in the TESS photometry. We apply
Bryant et al. (2020)’s method to our NGTS data assuming a circular Gaussian
PSF ($D$ = $K$ = 0) and a FWHM = 2 $\pm$ 0.25 pixels, so $S$ = 1.38 $\pm$ 0.3.
We determine a dilution factor of 0.64${}^{+0.33}_{-0.21}$% for our NGTS
photometry. Both of these dilution factors are used as inputs into the
EXOFASTv2 global model. We note that EXOFASTv2 defines their dilution factor
as the fractional contribution from neighboring stars, or $C$/(1+$C$).

##  4 Discussion

Our analysis finds NGTS-13b to have a mass of 4.84 $\pm$ 0.44 $\mathrm{M_{\rm
Jup}}$ and radius of 1.142 $\pm$ 0.046 $\mathrm{R_{\rm Jup}}$ at a short 4.12
day orbital period yielding an equilibrium temperature of 1605 $\pm$ 30 K,
making NGTS-13b a massive hot Jupiter. Although close-in very massive planets
are the easiest to detect, they are far from the most common. As of 16
September 2020 the NASA Exoplanet
Archive333http://exoplanetarchive.ipac.caltech.edu lists 4277 confirmed
planets; however, only 1764 have measured masses and only 1691 of these are
below 13 $\mathrm{M_{\rm Jup}}$ (the approximate border between planets and
brown dwarfs). Of planets with measured masses 971 are giant ($M_{P}$ = 0.5 -
13 $\mathrm{M_{\rm Jup}}$) with 223 of them having masses greater than 4
$\mathrm{M_{\rm Jup}}$. Here we only focus on planets that have well
determined masses and radii (percentage errors ¡ 10%), which consists of 225
planets, of which 155 are giants including only 17 massive giants ($M_{P}$ ¿ 4
$\mathrm{M_{\rm Jup}}$). These very massive planets have a similarly low
occurrence rate to the even more massive brown dwarfs, which are indeed known
to be rare (the brown dwarf desert, e.g., Marcy & Butler, 2000; Grether &
Lineweaver, 2006).

Figure 6: [Fe/H] and mass relationship for giant planets (0.5-13
$\mathrm{M_{\rm Jup}}$) with percentage errors $<$10% for both mass and radius
and a defined [Fe/H] in the NASA Exoplanet Archiverefnote1, which consists of
137 previously known planets including 17 with masses greater than 4
$\mathrm{M_{\rm Jup}}$. Figure 7: Radius (top) and density (bottom) as a
function of mass for transiting brown dwarfs and giant planets (0.5-13
$\mathrm{M_{\rm Jup}}$) with percentage errors $<$10% for both mass and radius
in the NASA Exoplanet Archiverefnote1. Transiting brown dwarfs are from the
literature summary found in (Mireles et al., 2020) for a total of 16 with
well-defined masses and radii. The stars display NGTS-13b and the vertical
dashed lines highlight 4 $\mathrm{M_{\rm Jup}}$ as well as 13 $\mathrm{M_{\rm
Jup}}$ (the approximate starting mass of brown dwarfs). The colors of the
symbols correspond to equilibrium temperature in the top plot and host star
mass in the bottom plot. The downward trend in maximum radii along increasing
mass exhibits the effect of stronger surface gravity in more massive objects,
which is clearly displayed by the linear trend in density. In the top plot we
only display planets that also have an equilibrium temperature in the Archive,
which consists of 133 previously known planets including 14 with masses
greater than 4 $\mathrm{M_{\rm Jup}}$.

###  4.1 Massive substellar companion populations

Previous studies have suggested that 4 $\mathrm{M_{\rm Jup}}$ could be an
approximate border for a distinct massive giant planet population compared to
lower-mass giant planets. Santos et al. (2017) find that massive giant planets
tend to orbit stars with $\left[{\rm Fe}/{\rm H}\right]$ distributions
statistically similar to field stars, unlike lower-mass giant planets which
occur more frequently around metal-rich stars. This could be evidence for two
different formation mechanisms where lower-mass giants may form via core
accretion (Pollack et al., 1996) and higher mass giants may form via disk
instability (Boss, 1997), as formation via disk instability is not as
metallicity dependent as core-accretion (Boss, 2002; Cai et al., 2006).

Additionally, Schlaufman (2018) find that substellar companions with masses
$\lesssim$ 4 $\mathrm{M_{\rm Jup}}$ preferentially orbit metal rich stars (a
property associated with core accretion), while companions with masses
$\gtrsim$ 10 $\mathrm{M_{\rm Jup}}$ do no share this property. Maldonado et
al. (2019) also find that massive giant planets tend to have host stars with
lower metallicities than lower-mass giant planets and conclude that the core-
accretion planet formation mechanism achieves its maximum efficiency for
planets with masses in the range 0.2 – 2 $\mathrm{M_{\rm Jup}}$. Goda & Matsuo
(2019) conclude that the mean metallicity of stars hosting companions between
4 and 25 $\mathrm{M_{\rm Jup}}$ is lower than that of the sample of companions
with masses 0.3 to 4 $\mathrm{M_{\rm Jup}}$, and the mean metallicity of stars
with companions more massive than 25 $\mathrm{M_{\rm Jup}}$ is much lower than
those of the other two sub-samples. Goda & Matsuo (2019) use 25
$\mathrm{M_{\rm Jup}}$ as the lower mass limit of brown dwarfs (instead of the
deuterium burning limit) because they found the mean metallicity of G-type
host stars with objects lighter than 25 $\mathrm{M_{\rm Jup}}$ to be much
higher than those with objects more massive than 25 $\mathrm{M_{\rm Jup}}$;
they suggest 25 $\mathrm{M_{\rm Jup}}$ is the upper mass limit for core-
accreted planets. However, other previous studies suggest lower mass brown
dwarfs (13 - 42.5 $\mathrm{M_{\rm Jup}}$) form via disk instability, e.g., Ma
& Ge (2014).

Adibekyan (2019) analyzes the known giant planet population and his results do
not support that 4 $\mathrm{M_{\rm Jup}}$ is a transition point between two
separate formation channels, but his results do suggest that high mass planets
can form through different mechanisms depending on their initial environment.
With a metallicity of $\left[{\rm Fe}/{\rm H}\right]$ = 0.25 $\pm$ 0.17
NGTS-13b’s host star is metal-rich and does not provide additional evidence
that more massive-giant planets prefer a different formation mechanism than
core accretion. We put NGTS-13b in this context in Figure 6 which shows the
metallicity of giant planets from 0.5 - 13 $\mathrm{M_{\rm Jup}}$.

###  4.2 NGTS-13b in context of the mass-radius relationship

We put NGTS-13b in terms of its density, $\rho_{P}$ = 4.02 $\pm$ 0.55 g cm-3,
in Figure 7 by displaying it in the mass-radius diagram with other giant
planets and brown dwarfs as well as density as a function of mass. We color
symbols in the mass-radius diagram based on equilibrium temperature to show
the effect high equilibrium temperatures can have on radii. We can see in
Figure 7 the larger maximum radii for planets with masses lower than $\sim$4
$\mathrm{M_{\rm Jup}}$, but this is explained by the redder coloring of these
symbols indicating higher equilibrium temperatures (resulting in more
inflation). NGTS-13b has a typical radius compared to planets with similar
masses suggesting that even with its high equilibrium temperature it is not
significantly affected by bloating (inflated radius due to a planet’s
atmosphere being heated from stellar irradiation and expanding), which is
expected for massive planets (Sestovic et al., 2018). The downward trend in
radii (and increasing trend in density) with larger mass in Figure 7
highlights the importance of the larger surface gravity in more massive
objects.

###  4.3 A slightly evolved host star

NGTS-13 has an effective temperature of $T_{\rm eff}$ = 5819 $\pm$ 73 K,
suggesting it is a G2 type star (Pecaut & Mamajek, 2013). However, the rather
large stellar radius of $R_{*}$ = 1.79 $\pm$ 0.06 $\mathrm{R_{\odot}}$
indicates NGTS-13 has evolved slightly as this radius is larger than expected
for a main sequence star of this temperature, and its surface gravity $\log
g_{*}$ = 4.04 $\pm$ 0.05 suggests the star is a subgiant ($\log g_{*}$ $<$
4.1). Our MIST isochrone analysis (see Section 3.2) finds an equivalent
evolutionary point (EEP) value of 413${}^{+47}_{-18}$, which still puts
NGTS-13 in the main sequence evolutionary phase, as the terminal age main
sequence EEP does not start until EEP = 454, see Choi et al. (2016), Dotter
(2016), and the MIST
documentation444http://waps.cfa.harvard.edu/MIST/README_tables.pdf.
Additionally, it is less likely in terms of evolutionary timescales to observe
a star near or past the turnoff point than to observe it in the middle of the
main sequence, therefore NGTS-13 is likely still fusing Hydrogen in its core
but has likely begun the transition from the main sequence to the red giant
branch placing it in the subgiant branch.

Lillo-Box et al. (2016) found a relative lack of hot Jupiters around giant and
subgiant stars, suggesting that the close-in massive planets around main-
sequence stars are engulfed by the star as it evolves. However, other studies
do not find a lack of hot Jupiters around evolved hosts; additionally NGTS-13
is likely still burning hydrogen in its core and has not evolved long enough
to engulf planets near the orbital distance of NGTS-13b. Grunblatt et al.
(2019) found that low-luminosity red giant branch evolved stars have a
tentatively higher population of close-in giant ($R_{P}$ $>$ 1 $\mathrm{R_{\rm
Jup}}$) planets than main-sequence stars suggesting that stellar evolution
does not affect close-in giant planet occurrence significantly until stars are
substantially in the red giant branch with radii of 5-6 $\mathrm{R_{\odot}}$.
Zhou et al. (2019) studied the occurrence rate of hot Jupiters as a function
of stellar mass and did not find any statistically consistent trends with
stellar mass and that hot Jupiters are just as abundant around main-sequence A
stars as they are around F and G stars.

The $M_{*}$ = 1.30${}^{+0.11}_{-0.18}$ $\mathrm{M_{\odot}}$ stellar mass of
NGTS-13 is larger than a typical G-type main-sequence star, which suggests it
evolved to a cooler temperature and is a ”retired” F star. A relatively larger
host star mass is consistent with previous studies that found massive planets
to be more common around massive stars (Santos et al., 2017; Maldonado et al.,
2019). Santos et al. (2017) find that more massive giant planets are more
common around more massive hosts, which are often more evolved; they suggest
the two distinct populations of giant planets could be related to evolved
stars not showing a clear metallicity-giant planet correlation.

We display host star mass by color in Figure 7. Only one host star has a mass
below 0.9 $\mathrm{M_{\odot}}$ in our sample of 17 well-defined massive giant
planets. Grunblatt et al. (2018) found that close-in giant planets around
evolved stars tend to have more eccentric orbits than those around main
sequence stars, but they focus on giant stars and orbital periods longer than
4.5 days. We explore NGTS-13b’s eccentricity further in Section 4.4.

###  4.4 A slightly eccentric orbit

Figure 8: Posterior distribution of the eccentricity from our EXOFASTv2 Markov
Chain Monte Carlo (MCMC) analysis including all chains and steps of the
distribution. Figure 9:  Eccentricity-period relationship for transiting brown
dwarfs and giant planets (0.5-13 $\mathrm{M_{\rm Jup}}$) with percentage
errors $<$10% for both mass and radius and a defined eccentricity in the NASA
Exoplanet Archiverefnote1, which consists of 129 previously known planets
including 17 with masses greater than 4 $\mathrm{M_{\rm Jup}}$. Transiting
brown dwarfs are from the literature summary found in (Mireles et al., 2020)
for a total of 16 with well-defined masses and radii.

NGTS-13b has a slightly eccentric, $e$ = 0.086 $\pm$ 0.034, orbit with a short
period of 4.12 days. Analyzing the posterior distribution of the eccentricity
displayed in Figure 8 we see that it is not centered around 0. We performed
the test of Lucy & Sweeney (1971), and find the statistical significance of
the eccentric fit to be P$(e>0)$ = 0.9591, which just passes the 5%
significance level suggested by Lucy & Sweeney (1971).

To further test the eccentric model we fit a Keplerian to the 13 RVs with
priors on the time of conjunction and period and compare fits with and without
forcing $\sqrt{e}\cos\omega_{*}$ and $\sqrt{e}\sin\omega_{*}$ to 0. We also
allow an RV jitter term to vary, which corresponds to $k$ = 7 free parameters
for the eccentric fit and $k$ = 5 for the circular fit. We determine the
Bayesian Information Criterion (BIC), where given the free parameters $k$,
number of measurements $n$, and maximized value of the likelihood function of
the model $\mathscr{L}$:

$\text{BIC}=k\ln n-2\ln\mathscr{L}.$ (2)

When picking from several models the one with the lowest BIC is preferred. We
find a $\Delta$BIC = 7.17 showing the eccentricity model is moderately favored
over the circularized model. For this test we found an RV jitter of 31 m s-1
for the circular fit but 0 m s-1 for the eccentric fit as no additional RV
uncertainty is needed for the eccentric model.

We note that for our final EXOFASTv2 analysis presented the model has a RV
jitter variance $\sigma_{J}^{2}$ = $-800^{+2400}_{-1100}$ which essentially
reduces the error size of the RVs when evaluating the model. As noted in
Equation 13 of Hara et al. (2019), the uncertainty in the eccentricity is
approximately proportional to the error of the RV measurement; therefore, this
added RV jitter variance may cause an underestimation of our eccentricity
uncertainties. Hara et al. (2019) show that the eccentricity estimate can be
affected by an undetected correlated signal. Given the amplitude of the RV
signal of NGTS-13b, we study the possibility that an undetected outer
planetary companion would cause a nonzero eccentricity. We test this by adding
a linear drift to our global EXOFASTv2 model. When allowing for an
unconstrained linear RV drift we find a slope of only 1.1${}^{+1.2}_{-1.4}$ m
s-1 day-1 and an eccentricity $e$ = 0.084${}^{+0.023}_{-0.021}$, similar to
our presented model with no RV drift. The similar eccentricity suggests that
NGTS-13b’s eccentricity measurement is not caused by the model attempting to
account for an undetected planet.

The slight eccentricity of NGTS-13b is not uncommon for massive planets (see
Figure 9). Hansen (2010) found that hot giant planets more massive than 3
$\mathrm{M_{\rm Jup}}$ have the upper envelope of their eccentricity
distribution shifted to lower periods compared to less massive hot Jupiters.
Also Dawson & Johnson (2018) found that most hot Jupiters with periods $<$ 3
days have circular orbits, but some hot-Jupiters in the 3-10 day orbital
period range occupy moderately eccentric orbits. Assuming a general tidal
circularization effect for close-in planets from equation 3 of Adams &
Laughlin (2006), NGTS-13b has a circularization timescale of 0.5 Gyr for a
tidal quality factor $Q_{p}$ = 105, 5.0 Gyr for $Q_{p}$ = 106, and 50.3 Gyr
for $Q_{p}$ = 107. Given the 4.23${}^{+2.65}_{-1.59}$ Gyr age of the system,
the current eccentricity may be an indication that NGTS-13b underwent
dynamical interactions with other components in the system during its
migration history, e.g., high-eccentricity tidal migration (see Section 3.1 of
Dawson & Johnson (2018)). We put NGTS-13b’s eccentricity into context with
other giant planets and transiting brown dwarfs in Figure 9.

##  5 Conclusions

We report the discovery of NGTS-13b, a massive hot Jupiter orbiting a subgiant
star at a 4.12 day period. NGTS-13b has a mass of 4.84 $\pm$ 0.44
$\mathrm{M_{\rm Jup}}$, radius of 1.142 $\pm$ 0.046 $\mathrm{R_{\rm Jup}}$ and
eccentricity $e$ = 0.086 $\pm$ 0.034 with an estimated equilibrium temperature
of 1605 $\pm$ 30 K. Even with this high equilibrium temperature NGTS-13b is a
dense object with $\rho_{P}$ = 4.02 $\pm$ 0.55 g cm-3 as is expected for
massive planets. The effective temperature, $T_{\rm eff}$ = 5819 $\pm$ 73 K,
of NGTS-13 suggests it is a G-type star; however, its mass of $M_{*}$ =
1.30${}^{+0.11}_{-0.18}$ $\mathrm{M_{\odot}}$ and radius $R_{*}$ = 1.788 $\pm$
0.057 $\mathrm{R_{\odot}}$ suggest it is a ”retired” F star in the subgiant
branch.

NGTS-13b is a valuable addition to the relatively rare group ($<$20 currently
known) of well-characterized (mass and radius percent errors $<$ 10%) massive
($M_{P}$ $>$ 4 $\mathrm{M_{\rm Jup}}$) giant planets that are needed to
understand the differences between giants, massive giants, and brown dwarfs.
Previous studies have considered 4 $\mathrm{M_{\rm Jup}}$ as a boundary
between core accretion planet formation and other formation mechanisms such as
disk instability, based on metallicity distributions. However, NGTS-13 has a
metallicity $\left[{\rm Fe}/{\rm H}\right]$ = 0.25 $\pm$ 0.17 and does not
give further evidence that these massive giants form through other mechanisms
than core-accretion.

###### Acknowledgements.

We thank the Swiss National Science Foundation (SNSF) and the Geneva
University for their continuous support to our planet search programs. This
work was carried out in the frame of the National Centre for Competence in
Research $PlanetS$ supported by the Swiss National Science Foundation (SNSF).
This work uses data collected under the NGTS project at the ESO La Silla
Paranal Observatory. The NGTS facility is operated by the consortium
institutes with support from the UK Science and Technology Facilities Council
(STFC) under projects ST/M001962/1 and ST/S002642/1. This publication makes
use of The Data & Analysis Center for Exoplanets (DACE), which is a facility
based at the University of Geneva (CH) dedicated to extrasolar planets data
visualisation, exchange and analysis. DACE is a platform of the Swiss National
Centre of Competence in Research (NCCR) PlanetS, federating the Swiss
expertise in Exoplanet research. The DACE platform is available at
https://dace.unige.ch. This paper includes data collected by the TESS mission.
Funding for the TESS mission is provided by the NASA Explorer Program. This
work has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC, https://www.cosmos.
esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by
national institutions, in particular the institutions participating in the
Gaia Multilateral Agreement. JSJ acknowledges support by FONDECYT grant
1201371 and partial support from CONICYT project Basal AFB-170002. JVS
acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (project Four
Aces; grant agreement No. 724427). MNG acknowledges support from MIT’s Kavli
Institute as a Juan Carlos Torres Fellow. Contributions by authors from the
University of Warwick were supported by STFC consolidated grants ST/P000495/1
and ST/T000406/1. EG gratefully acknowledges support from the David and
Claudia Harding Foundation in the form of a Winton Exoplanet Fellowship. The
authors thank Aaron Dotter for his helpful advice.

## References

  * Adams & Laughlin (2006) Adams, F. C. & Laughlin, G. 2006, ApJ, 649, 1004
  * Adibekyan (2019) Adibekyan, V. 2019, Geosciences, 9, 105
  * Alam et al. (2015) Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12
  * Barker & Ogilvie (2009) Barker, A. J. & Ogilvie, G. I. 2009, MNRAS, 395, 2268
  * Bianchi et al. (2011) Bianchi, L., Herald, J., Efremova, B., et al. 2011, Ap&SS, 335, 161
  * Boss (1996) Boss, A. P. 1996, in Lunar and Planetary Science Conference, Vol. 27, Lunar and Planetary Science Conference, 139
  * Boss (1997) Boss, A. P. 1997, Science, 276, 1836
  * Boss (2002) Boss, A. P. 2002, ApJ, 567, L149
  * Bouma et al. (2019) Bouma, L. G., Hartman, J. D., Bhatti, W., Winn, J. N., & Bakos, G. Á. 2019, ApJS, 245, 13
  * Bryant et al. (2020) Bryant, E. M., Bayliss, D., Nielsen, L. D., et al. 2020, arXiv e-prints, arXiv:2009.10620
  * Cai et al. (2006) Cai, K., Durisen, R. H., Michael, S., et al. 2006, ApJ, 636, L149
  * Carmichael et al. (2020a) Carmichael, T. W., Quinn, S. N., Mustill, A. J., et al. 2020a, AJ, 160, 53
  * Carmichael et al. (2020b) Carmichael, T. W., Quinn, S. N., Zhou, G., et al. 2020b, arXiv e-prints, arXiv:2009.13515
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Claret & Bloemen (2011) Claret, A. & Bloemen, S. 2011, A&A, 529, A75
  * Cubillos et al. (2017) Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3
  * D’Angelo et al. (2003) D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540
  * Dawson & Johnson (2018) Dawson, R. I. & Johnson, J. A. 2018, ARA&A, 56, 175
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * Eastman (2017) Eastman, J. 2017, EXOFASTv2: Generalized publication-quality exoplanet modeling code
  * Eastman et al. (2013) Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83
  * Eastman et al. (2019) Eastman, J. D., Rodriguez, J. E., Agol, E., et al. 2019, arXiv e-prints, arXiv:1907.09480
  * Ford (2006) Ford, E. B. 2006, ApJ, 642, 505
  * Gelman et al. (2003) Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 2003, Bayesian Data Analysis, 2nd edn. (London: Chapman & Hall)
  * Gelman & Rubin (1992) Gelman, A. & Rubin, D. B. 1992, Statistical Science, 7, 457
  * Gill et al. (2019) Gill, S., Maxted, P. F. L., Evans, J. A., et al. 2019, A&A, 626, A119
  * Gill et al. (2018) Gill, S., Maxted, P. F. L., & Smalley, B. 2018, A&A, 612, A111
  * Gill et al. (2020) Gill, S., Wheatley, P. J., Cooke, B. F., et al. 2020, ApJ, 898, L11
  * Goda & Matsuo (2019) Goda, S. & Matsuo, T. 2019, ApJ, 876, 23
  * Grether & Lineweaver (2006) Grether, D. & Lineweaver, C. H. 2006, ApJ, 640, 1051
  * Grunblatt et al. (2019) Grunblatt, S. K., Huber, D., Gaidos, E., et al. 2019, AJ, 158, 227
  * Grunblatt et al. (2018) Grunblatt, S. K., Huber, D., Gaidos, E., et al. 2018, ApJ, 861, L5
  * Hansen (2010) Hansen, B. M. S. 2010, ApJ, 723, 285
  * Hara et al. (2019) Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738
  * Henden et al. (2016) Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog, II/336
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Ikwut-Ukwa et al. (2020) Ikwut-Ukwa, M., Rodriguez, J. E., Bieryla, A., et al. 2020, AJ, 160, 209
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  * Lee & Chiang (2016) Lee, E. J. & Chiang, E. 2016, ApJ, 817, 90
  * Lee et al. (2014) Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95
  * Lillo-Box et al. (2016) Lillo-Box, J., Barrado, D., & Correia, A. C. M. 2016, A&A, 589, A124
  * Lin et al. (1996) Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606
  * Lucy & Sweeney (1971) Lucy, L. B. & Sweeney, M. A. 1971, AJ, 76, 544
  * Ma & Ge (2014) Ma, B. & Ge, J. 2014, MNRAS, 439, 2781
  * Maldonado et al. (2019) Maldonado, J., Villaver, E., Eiroa, C., & Micela, G. 2019, A&A, 624, A94
  * Marcy & Butler (2000) Marcy, G. W. & Butler, R. P. 2000, PASP, 112, 137
  * Mayor & Queloz (1995) Mayor, M. & Queloz, D. 1995, Nature, 378, 355
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
  * Mireles et al. (2020) Mireles, I., Shporer, A., Grieves, N., et al. 2020, arXiv e-prints, arXiv:2006.14019
  * Pál (2009) Pál, A. 2009, PhD thesis, Department of Astronomy, Eötvös Loránd University
  * Pecaut & Mamajek (2013) Pecaut, M. J. & Mamajek, E. E. 2013, ApJS, 208, 9
  * Pepe et al. (2002) Pepe, F., Mayor, M., Rupprecht, G., et al. 2002, The Messenger, 110, 9
  * Pepper et al. (2020) Pepper, J., Kane, S. R., Rodriguez, J. E., et al. 2020, AJ, 159, 243
  * Piso et al. (2015) Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Queloz et al. (2001) Queloz, D., Mayor, M., Udry, S., et al. 2001, The Messenger, 105, 1
  * Rafikov (2006) Rafikov, R. R. 2006, ApJ, 648, 666
  * Rasio & Ford (1996) Rasio, F. A. & Ford, E. B. 1996, Science, 274, 954
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
  * Santos et al. (2017) Santos, N. C., Adibekyan, V., Figueira, P., et al. 2017, A&A, 603, A30
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schlaufman (2018) Schlaufman, K. C. 2018, ApJ, 853, 37
  * Sestovic et al. (2018) Sestovic, M., Demory, B.-O., & Queloz, D. 2018, A&A, 616, A76
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Stassun et al. (2019) Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138
  * Stassun & Torres (2016) Stassun, K. G. & Torres, G. 2016, AJ, 152, 180
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466
  * Triaud et al. (2010) Triaud, A. H. M. J., Collier Cameron, A., Queloz, D., et al. 2010, A&A, 524, A25
  * Wheatley et al. (2018) Wheatley, P. J., West, R. G., Goad, M. R., et al. 2018, MNRAS, 475, 4476
  * Winn et al. (2005) Winn, J. N., Noyes, R. W., Holman, M. J., et al. 2005, ApJ, 631, 1215
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Yee et al. (2017) Yee, S. W., Petigura, E. A., & von Braun, K. 2017, ApJ, 836, 77
  * Zhou et al. (2019) Zhou, G., Huang, C. X., Bakos, G. Á., et al. 2019, AJ, 158, 141


# Robust Tracking Control for Nonlinear Systems: Performance optimization via
extremum seeking

Jiapeng Xu1, Ying Tan2, and Xiang Chen1 1Jiapeng Xu and Xiang Chen are with
the Department of Electrical and Computer Engineering, University of Windsor,
Windsor, ON N9B 3P4, Canada<EMAIL_ADDRESS><EMAIL_ADDRESS>Tan is
with the Department of Mechanical Engineering, The University of Melbourne,
Melbourne, VIC 3010, Australia<EMAIL_ADDRESS>
###### Abstract

This paper presents a controller design and optimization framework for
nonlinear dynamic systems to track a given reference signal in the presence of
disturbances when the task is repeated over a finite-time interval. This novel
framework mainly consists of two steps. The first step is to design a robust
linear quadratic tracking controller based on the existing control structure
with a Youla-type filter $\tilde{Q}$. Secondly, an extra degree of freedom: a
parameterization in terms of $\tilde{Q}$, is added to this design framework.
This extra design parameter is tuned iteratively from measured tracking cost
function with the given disturbances and modeling uncertainties to achieve the
best transient performance. The proposed method is validated with simulation
placed on a Furuta inverted pendulum, showing significant tracking performance
improvement.

##  I Introduction

Robust tracking control of nonlinear systems has been extensively studied in
the literature using various robust techniques, such as $H_{\infty}$ control
[1, 2] and sliding mode control [3]. These methods are in general the worst-
case design, which would ensure the stability under the worst-case
disturbances. On the other hand, optimal performance such as a linear
quadratic form has been the focus of the design in industry applications.
However, it is usually hard to analyze the performance for nonlinear dynamics.
One of the key reasons to contribute this difficulty comes from the the
analysis tool used in stability analysis. Lyapunov direct method [4] has been
used to provide sufficient conditions to guarantee the stability of nonlinear
dynamics. The optimal control for nonlinear dynamics requires to solve the
Hamilton–Jacobi–Bellman (HJB) equation, which is a nonlinear partial
differential equation with respect to a given cost. Solving this HJB equation
is computational costly.

In contrast, both robust and optimal control designs have been extensively
investigated for linear time-invariant (LTI) dynamic systems [5, 6, 7, 8]. In
particular, a robust controller design with a Youla-type filter $\tilde{Q}$
[8, 9] has been proposed recently, which is motivated by the generalized
internal model control (GIMC) proposed in [6]. The robust controller with
$\tilde{Q}$ provides automatic robustness recovery in the linear quadratic
Gaussian (LQG)/$H_{2}$ control [8]. Its key idea is to use the $\tilde{Q}$
filter to balance the optimal performance without consideration of the
disturbance and robust performance using the techniques such as $H_{\infty}$
control. Since the filter $\tilde{Q}$ is driven by the residual signal
indicating the mismatch between the nominal model and the true system,
$\tilde{Q}$ is only activated when there exists unmodelled dynamics or
external disturbances, such that this kind of controller design can lead to a
high performance in the presence of disturbances and uncertainties. This
technique is quite different from the traditional mixed $H_{2}/H_{\infty}$
control, a trade-off design [10, 11, 12].

This work proposes to utilizes the robust controller with $\tilde{Q}$ in [8]
to systematically design the feedback control for a nonlinear dynamic system
via its linearization. The proposed framework is used to generate optimal
tracking performance for a class of nonlinear dynamic systems to track a given
reference trajectory. More specifically, the proposed framework first presents
a robust linear quadratic tracking (LQT) controller design based on the filter
$\tilde{Q}$. Then by introducing an extra gain factor in terms of $\tilde{Q}$,
which can be treated as the balance between the LQT performance and the
robustness with respect to disturbances and uncertainties coming from
linearizations and other external signals, an updating law is generated to
tune this gain factor to minimize the tracking cost in the presence of
modeling uncertainties and disturbances. The choice of the gain factor does
not affect the local stability properties of the closed-loop nonlinear system,
while it improves the tracking performance. In this work, the data-driven
extremum seeking (ES) approach [13, 14, 15, 16], which is a model-free
optimization method, is adapted to find this optimal gain factor.
Alternatively, other model-free optimization techniques such as reinforcement
learning can also be considered [2, 17].

The effectiveness of the proposed framework is validated with simulation
placed on a Furuta inverted pendulum. It has been shown that this optimal
parameter is dependent on the nonlinear dynamics, the type of the reference
trajectories, as well as the type of disturbances. The obtained optimal gain
can achieve much better transient tracking performance compared with the
standard LQT controller and the robust controllers such as $H_{\infty}$.

The remainder of this paper is organized as follows. Section II formulates the
tracking problem for nonlinear systems of interest. Section III presents a
design procedure of LQT controller with $\tilde{Q}$. Section IV further
presents a controller design for nonlinear systems where performance is
further optimized via ES. Section V provides simulation results on an inverted
pendulum. Finally, Section VI concludes this work.

##  II Problem Formulation

Nonlinear systems of the following form are considered:

$\displaystyle\dot{x}$ $\displaystyle=f(x,u,w),\hskip
14.45377ptx(0)\in{\mathcal{R}}^{n}$ $\displaystyle y$ $\displaystyle=g(x,w),$
(1)

where $x\in{\mathcal{R}}^{n}$ is the state of system, $y\in{\mathcal{R}}^{p}$
is the measurement output, $u\in{\mathcal{R}}^{m}$ is the control input, and
$w\in{\mathcal{R}}^{n_{w}}$ is a disturbance signal representing external
disturbances and/or modeling uncertainties. The nonlinear mappings
$f:{\mathcal{D}}_{X}\times{\mathcal{D}}_{U}\times{\mathcal{D}}_{W}\rightarrow{\mathcal{R}}^{n}$
and $g:{\mathcal{R}}^{n}\times{\mathcal{D}}_{W}\rightarrow{\mathcal{R}}^{p}$
are continuously differentiable functions. Here ${\mathcal{D}}_{X}$,
${\mathcal{D}}_{U}$, and ${\mathcal{D}}_{W}$ are compact sets in
${\mathcal{R}}^{n}$, ${\mathcal{R}}^{m}$, and ${\mathcal{R}}^{n_{w}}$
respectively. Moreover, it is assumed that $f(0,0,0)=0$, indicating the origin
is in the set
$\Omega={\mathcal{D}}_{X}\times{\mathcal{D}}_{U}\times{\mathcal{D}}_{W}$.

The linearization of system (1) around $(x=0,u=0,w=0)$ becomes the following
LTI system:

$\displaystyle\dot{x}$ $\displaystyle=Ax+B_{1}w+B_{2}u,$ (2) $\displaystyle y$
$\displaystyle=C_{2}x+D_{21}w,$ (3)

where

$\displaystyle A$ $\displaystyle=\frac{\partial f}{\partial
x}\bigg{|}_{x,u,w=0},\;B_{1}=\frac{\partial f}{\partial
w}\bigg{|}_{x,u,w=0},\;\;B_{2}=\frac{\partial f}{\partial
u}\bigg{|}_{x,u,w=0},$ $\displaystyle C_{2}$ $\displaystyle=\frac{\partial
g}{\partial x}\bigg{|}_{x,w=0},\;\;D_{21}=\frac{\partial g}{\partial
w}\bigg{|}_{x,w=0}.$

The following assumption is standard to stabilize the linearized system (2) by
output feedback using the output signal in (3).

###### Assumption 1

(i) $(A,B_{2})$ is stabilizable, and (ii) $(C_{2},A)$ is detectable.

The control objective of this work is to design an appropriate controller such
that the output

$\displaystyle\tilde{y}=Ex,\;E\in{\mathcal{R}}^{p_{1}\times n},$ (4)

can track a desired (reference) trajectory $r$ that is known _a priori_.
Specifically, we seek to find an optimal control input to minimize the
following finite-time quadratic cost function:

$\displaystyle
J(u)=\frac{1}{T}\int_{0}^{T}[(\tilde{y}-r)^{\prime}Q(\tilde{y}-r)+u^{\prime}Ru]{\rm
d}t,$ (5)

where $Q$ is positive semi-definite and $R$ is positive definite. In this
paper, we use “′” to represent the transpose of a matrix.

We propose to solve the formulated tracking problem for nonlinear systems (1)
by two steps. The first step is to design a robust linear quadratic tracking
(LQT) controller based on a Youla-type filter $\tilde{Q}$ proposed in [8] ,
which will be detailed in the next section. The second step is to introduce an
extra gain factor to the filter $\tilde{Q}$ and use the data-driven
optimization algorithm such as extremum seeking (ES) [13, 14, 15, 16] to tune
this factor.

##  III Robust Linear Quadratic Tracking Control with $\tilde{Q}$

This section presents a robust LQT controller design based on the filter
$\tilde{Q}$ in [8], which balances the robustness with respect to disturbances
and the LQ optimal tracking control, for the LTI systems taking the same form
as the linearized system (2) and (3).

###  III-A Preliminaries of the controller structure with $\tilde{Q}$

The structure of the controller with $\tilde{Q}$ is presented in Fig. 1 (see
more discussions in [8]). It contains some performance related variable $z$
defined as

$\displaystyle z=C_{1}x+D_{12}u,$ (6)

an observer, and the filter $\tilde{Q}(\cdot)$ satisfying $\tilde{Q}(0)=0$.
The input of the filter $\tilde{Q}$ is the residual signal $f=\hat{y}-y$,
which is the deviation of estimated and actual sensor outputs. It reflects the
mismatch between the nominal model and the true system. As shown in [8], the
output of the filter $\tilde{Q}(\cdot)$ is designed for robustness recovery,
reducing the effect of disturbances. For example, the $H_{\infty}$ design
method can be used.

The control input $u$ consists of two parts : $u_{l}$ and $u_{f}$. The first
part $u_{l}$ is a nominal control, which is to design the stabilizing
controller with the help of the observer without consideration of $w(t)$. The
matrices $(F,L)$ are the state feedback and observer gains, respectively, such
that $A+B_{2}F$ and $A+LC_{2}$ are Hurwitz by Assumption 1. The second part
$u_{f}$ is the output of the filter, i.e., $u_{f}=\tilde{Q}(f)$. If
$w(t)=0,\forall t\geq 0$, the residual signal $f$ is zero and produces
$u_{f}=\tilde{Q}(f)=0$. Hence the role of this filter is to enhance the
robustness. One admissible filter $\tilde{Q}$ is of the following form [8]:

$\displaystyle\tilde{Q}:\;\dot{x}_{q}=A_{q}x_{q}+B_{q}f,\;\;u_{f}=F_{q}x_{q}.$
(7)

Matrices $A_{q}$, $B_{q}$ and $F_{q}$ can be obtained by the standard
$H_{\infty}$ techniques [18, 5] for the augmented system with state $(x,e)$,
with $e=x-\hat{x}$, and $A_{q}$ is Hurwitz.

Figure 1: The diagram of the controller with $\tilde{Q}$.

Following this controller structure, next will introduce the LQT controller
design, followed by the robust controller using $\tilde{Q}$.

###  III-B Linear Quadratic Tracking (LQT)

Consider the nominal case of the LTI system (2) and (3), i.e., $w=0$ and the
initial state $x(0)$ is known:

$\displaystyle\dot{x}$ $\displaystyle=Ax+B_{2}u,\;x(0)~{}\text{given}.$ (8)

This reduces to a state feedback control problem. The following proposition
provides the optimal solution to minimizing (5) by using the Hamilton-Jacobi-
Bellman equation [7, Chapter 2]; see also [7, Chapter 4].

###### Proposition 1

Considering the system consisting of (8) and (4), for a given reference
trajectory $r$ and the given performance index (5), the optimal control
$u^{*}$ is given by

$\displaystyle-\dot{b}$
$\displaystyle=(A+B_{2}F)^{\prime}b+E^{\prime}Qr,\;b(T)=0,$ (9) $\displaystyle
u^{*}$ $\displaystyle={Fx}+{R^{-1}B_{2}^{\prime}b},$ (10)

where $F=-R^{-1}B_{2}^{\prime}P$, and $P$ is the solution of the following
differential Riccati equation:

$\displaystyle-\dot{P}=PA+A^{\prime}P-$ $\displaystyle
PB_{2}R^{-1}B_{2}^{\prime}P+E^{\prime}QE.$ (11)

In practice, to simplify the design of LQT controller, the algebraic Riccati
equation instead of differential Riccati equation (11) can be considered, such
that a time-invariant stabilizing control gain $F$ is obtained. For
convenience, the feed-forward control law in (10) is denoted as
$u_{r}:=R^{-1}B_{2}^{\prime}b$. The diagram of control system with LQ optimal
tracking is shown in Fig. 2.

Figure 2: The diagram of LQ optimal tracking control.

###  III-C Design of LQT with $\tilde{Q}$

Now we design the LQT controller with $\tilde{Q}$ to robustify the feedback
system. To do this, we redefine the performance variable $z$ in (6) as

$\displaystyle z=C_{1}(\tilde{y}-r)+D_{12}u,$ (12)

such that it coincides with the tracking cost in (5). We can let
$C_{1}^{\prime}C_{1}=Q$, $D_{12}^{\prime}D_{12}=R$ and
$C_{1}^{\prime}D_{12}=0$ to reduce the weighing matrix parameters. To design
the filter $\tilde{Q}$, set the reference signal $r$ to be zero such that
$u_{r}=0$ and introduce an observer gain $L$ such that $A+LC_{2}$ is Hurwitz.
Then the design procedure of $\tilde{Q}$ can be the same as in [8], i.e., (7).
More specifically, the filter $\tilde{Q}$ is designed according to the
following $H_{\infty}$ performance criterion:

$\displaystyle\|T_{zw}(s)\|_{\infty}<\gamma,$ (13)

where $\gamma>0$ is a prescribed value and $T_{zw}(s)$ is the closed-loop
transfer function from $w$ to $z$. See [8] for more details.

By combing the LQT controller, the observer gain $L$ and the filter
$\tilde{Q}$ designed using $H_{\infty}$ technique, the diagram of the robust
LQT with $\tilde{Q}$ is shown in Fig. 3.

Figure 3: The diagram of robust LQT with parameterized $\tilde{Q}$.

###### Remark 1

The transfer functions from $u_{r}$ to $u$ ,$y$, and $\tilde{y}$ in the
closed-loop system are independent of the filter $\tilde{Q}$. These transfer
function are the same as those in the LQT design. Therefore, the proposed
tracking design scheme allows us to select $\tilde{Q}$ to attenuate
disturbances without affecting the tracking capability in the LQT design.

##  IV Performance Optimization for Nonlinear Systems via Extremum Seeking

We propose to solve the formulated optimal tracking problem based on the
linearized model (2) and (3). However, the performance of the robust LQT
controller with $\tilde{Q}$ designed in the previous section cannot be
guaranteed for nonlinear systems due to the existence of linearization errors.
Moreover, this linear controller cannot give an optimal control input
minimizing the performance index (5) for our tracking problem for nonlinear
systems (1). In the sequel, we will design an optimal control input by
introducing an extra gain factor $\alpha\in{\mathcal{R}}$ to the filter
$\tilde{Q}$ and then optimize this factor by the data-driven ES algorithm.

###  IV-A An extra gain factor for $\tilde{Q}$

To deal with linearization errors and unmodeled disturbances, a gain factor
$\alpha\in{\mathcal{R}}$ for the output of the filter $\tilde{Q}$ is
introduced as an extra degree of freedom, such that a new $\tilde{Q}$ denoted
by $\tilde{Q}_{\alpha}$ is obtained:

$\displaystyle\tilde{Q}_{\alpha}:\;\dot{x}_{q}=A_{q}x_{q}+B_{q}f,\;\;u_{f}=\alpha
F_{q}x_{q}.$ (14)

When $\alpha=0$, it leads to an optimal tracking controller as indicated in
Proposition 1 while $\alpha=1$ leads to the LQT controller with $\tilde{Q}$
shown in Fig. 3. Intuitively, the choice of this $\alpha$ reflects some
balance between the LQT and the robustness. This leads to a closed-loop
nonlinear system with the LQT controller with $\tilde{Q}_{\alpha}$ shown in
Fig. 4.

Now the tracking performance index $J(u)$ in (5) can be written as a function
of $\alpha$, i.e.,

$\displaystyle
J(\alpha)=\frac{1}{T}\int_{0}^{T}[(\tilde{y}-r)^{\prime}Q(\tilde{y}-r)+u^{\prime}Ru]{\rm
d}t,$ (15)

and we seek to find an optimal $\alpha$ parameter to minimize the above cost
function. Before solving the $\alpha$ optimization problem, we present the
following result about the local stability properties of such a closed-loop
nonlinear system.

Figure 4: The diagram of robust LQT control with parameterized
$\tilde{Q}_{\alpha}$ for nonlinear systems.

###### Theorem 1

Let the origin ($x=0,u=0,w=0$) be an equilibrium point for the nonlinear
system (1) and let system (2) and (3) be the linearization of (1) about the
origin. Then for any $\alpha\in{\mathcal{R}}$, the origin of the closed-loop
nonlinear system in Fig. 4 is locally exponentially stable.

###### Proof:

Let
$x_{c}:=\left[\begin{array}[]{cc}x_{q}^{\prime}&\hat{x}^{\prime}\end{array}\right]^{\prime}$.
Then it follows from the design of the LQT controller with $\tilde{Q}$ shown
in Fig. 3 and the expression of $\tilde{Q}_{\alpha}$ in (14) that we have the
following linear tracking controller:

$\displaystyle\dot{x}_{c}$
$\displaystyle=A_{c}x_{c}+B_{c}y+B_{r}u_{r},\;u=F_{c}x_{c}+u_{r},$ (16)

where

$\displaystyle A_{c}$
$\displaystyle=\left[\begin{array}[]{cc}A_{q}&B_{q}C_{2}\\\ \alpha
B_{2}F_{q}&A+B_{2}F+LC_{2}\end{array}\right],\;B_{c}=\left[\begin{array}[]{c}-B_{q}\\\
-L\end{array}\right],$ $\displaystyle B_{r}$
$\displaystyle=\left[\begin{array}[]{c}0\\\
B_{2}\end{array}\right],\;F_{c}=\left[\begin{array}[]{cc}\alpha
F_{q}&F\end{array}\right].$

Hence, by a linear transformation, it can be shown that the closed-loop matrix
for the linearized system (2) and (3) has the following triangular form

$\displaystyle\bar{A}=\bar{A}(\alpha)=\left[\begin{array}[]{ccc}A+B_{2}F&\alpha
B_{2}F_{q}&B_{2}F\\\ 0&A_{q}&B_{q}C_{2}\\\ 0&0&A+LC_{2}\end{array}\right].$
(20)

Since $A+B_{2}F$, $A+LC_{2}$ and $A_{q}$ are all Hurwitz, $\bar{A}$ is
Hurwitz. It can be verified that the origin $(x=0,x_{c}=0,w=0)$ is an
equilibrium point of the closed-loop nonlinear system. Thus, we can conclude
the local exponential stability of the closed-loop nonlinear system in Fig. 4
at the origin for any $\alpha\in{\mathcal{R}}$[4, Section 12.2].

Theorem 1 shows that the choice of $\alpha$ does not affect the local
stability properties of the closed-loop nonlinear system presented in Fig. 4
while it affects the performance, in particular, the transient performance.
This work tries to find an optimal $\alpha$ to tracking control for nonlinear
systems. In the sequel, the data-driven extremum seeking (ES) [13, 14, 15, 16]
will be used tune the parameter $\alpha$ for a given disturbance $w$.

###  IV-B Performance optimization via extremum seeking

It is highlighted that in the linearized model consisting of (2) and (3), the
disturbances coming from two parts: one is unmodelled uncertainties coming
from the linearization residue and unmodelled dynamics. Such uncertainties are
related to the size of compact sets ${\mathcal{D}}_{X}$, ${\mathcal{D}}_{U}$,
and ${\mathcal{D}}_{W}$. The other is other types of deterministic and random
noises. In this paper, our focus is the unmodelled uncertainties and
deterministic noises that are repeatable when the nonlinear dynamics run over
a fixed time interval $[0,T]$. Thus, the disturbance can be re-written as
$w=w(t,x,u)$.

###### Assumption 2

It is assumed that for a deterministic and repeatable disturbance
$w(t,x,u),t\in[0,T]$, for any $x\in{\mathcal{D}}_{X}$, and
$u\in{\mathcal{D}}_{U}$, there exists a unique optimal
$\alpha^{*}\in{\mathcal{R}}$ such that the cost function defined (15) can
reach a minimum.

Under this assumption, the ES algorithm presented in [13] is used to tune the
parameter $\alpha$ for a given disturbance $w(t,x,u)$ by repeatedly running
the closed-loop nonlinear system consisting of (1) with the control law (16)
over the finite time $[0,T]$. ES is a model-free optimization method which
uses only input–output data to see an optimal input with respect to a given
cost [13]. The ES algorithm adopted here works for the iteration domain, which
is the same as the one used in [14]:

$\displaystyle\zeta(k+1)$ $\displaystyle=-h\zeta(k)+J(\alpha(k)),$
$\displaystyle\hat{\alpha}(k+1)$
$\displaystyle=\hat{\alpha}(k)-\delta\beta\cos(\omega
k)[J(\alpha(k))-(1+h)\zeta(k)],$ $\displaystyle\alpha(k+1)$
$\displaystyle=\hat{\alpha}(k+1)+\beta\cos(\omega(k+1)),$ (21)

where $k$ is the iteration number, $0<h<1$, $\zeta(k)$ is a scalar, $\delta$
is the step size, and $\beta$ is the perturbation amplitude. Stability and
convergence are mainly influenced by the values of $\delta$ and $\beta$. The
modulation frequency $\omega$ is chosen such that $\omega=a\pi$, where $a$
satisfies $0<a<1$. The overall ES $\alpha$ tuning scheme is summarized in Fig.
5. The local convergence analysis of such ES algorithm when Assumption 2 holds
locally can be found in [19]. Using the similar analysis techniques as in [15]
for the continuous-time systems, non-local convergence can be achieved when
Assumption 2 holds.

The design procedure of the proposed robust optimal tracking control scheme
can be summarized as the following three steps.

  * •

Step 1: Obtain the linearized model around the origin $(x=0,u=0,w=0)$ from
nonlinear system (1);

  * •

Step 2: Design the LQT controller with $\tilde{Q}_{\alpha}$ as in Fig. 3 in
which $\tilde{Q}$ in (7) is replaced by $\tilde{Q}_{\alpha}$ in (14);

  * •

Step 3: Tune $\alpha$ via the ES algorithm (21) for a given disturbance
$w(t,x,u)$, shown in Fig. 5.

Figure 5: The overall ES $\alpha$ tuning scheme. The ES algorithm updates the
parameter $\alpha(k)$ in the filter $\tilde{Q}_{\alpha}$ to minimize the cost
$J(\alpha)$.

###### Remark 2

It is noted that the focus of this work is to ensure that the transient
behaviour of the closed-loop system presented in Fig. 5 over the finite time
interval $[0,T]$ can be improved over iteration when the disturbance $w$ is
repeatable over iteration. If Assumption 2 holds, the ES algorithm in (21) can
ensure that the cost function (15) decreases over iterations for appropriately
tuned parameters $(a,\delta,\beta,h)$ in (21). The decrease of the cost
function has been verified in simulations.

##  V Simulation Results

In order to verify the effectiveness of the proposed robust tracking control
in Fig. 5, a simulation example on a Furula inverted pendulum (the Quanser
QUBE Servo 2 rotary pendulum [20]) is presented.

###  V-A Linearization and design parameters

The Furuta inverted pendulum is an under-actuated (unstable) nonlinear system
with the state variable
$x=\left[\begin{array}[]{cccc}\theta_{1}&\theta_{2}&\dot{\theta}_{1}&\dot{\theta}_{2}\end{array}\right]^{\prime}$,
consisting of the base arm angle $\theta_{1}$ and the vertical arm (pendulum)
angle $\theta_{2}$, as well as the corresponding angular velocities. When the
pendulum is in the upright position,
$x=\left[\begin{array}[]{cccc}0&0&0&0\end{array}\right]^{\prime}$. The
nonlinear equations of motion for the pendulum system in disturbance-free case
are [20]:

$\displaystyle M\ddot{\theta}_{1}=$ $\displaystyle
m_{p}lr\ddot{\theta}_{2}\cos\theta_{2}-J_{p}\dot{\theta}_{1}\dot{\theta}_{2}\sin
2\theta_{2}-m_{p}lr\dot{\theta}_{2}^{2}\sin\theta_{2}$ $\displaystyle-
b_{r}\dot{\theta}_{1}+\tau,$ $\displaystyle J_{p}\ddot{\theta}_{2}=$
$\displaystyle
m_{p}lr\ddot{\theta}_{1}\cos\theta_{2}+0.5J_{p}\dot{\theta}_{1}^{2}\sin
2\theta_{2}+m_{p}gl\sin\theta_{2}$ $\displaystyle-b_{p}\dot{\theta}_{2},$ (22)

where $M=J_{r}+J_{p}\sin^{2}\theta_{2}$ and
$\tau=k_{m}/R_{m}(u-k_{m}\dot{\theta}_{1})$ is the applied torque at the base
of the rotary arm with $u$ the control input/voltage.

Following the 3-step design procedure in the previous section, we need to
linearize the nonlinear system first. The origin ($x=0,u=0$) is an equilibrium
point for pendulum system (V-A), leading to the following linearization model
with assumed disturbance/noise $w$

$\displaystyle A$ $\displaystyle=\left[\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&1\\\
0&\frac{m_{p}^{2}l^{2}rg}{J_{t}}&-\frac{J_{p}b_{r}}{J_{t}}-\frac{k_{m}^{2}J_{p}}{R_{m}J_{t}}&-\frac{m_{p}lrb_{p}}{J_{t}}\\\
0&\frac{J_{r}m_{p}gl}{J_{t}}&-\frac{m_{p}lrb_{r}}{J_{t}}-\frac{k_{m}^{2}m_{p}lr}{R_{m}J_{t}}&-\frac{J_{r}b_{p}}{J_{t}}\end{array}\right],$
$\displaystyle B_{2}$ $\displaystyle=\left[\begin{array}[]{c}0\\\ 0\\\
\frac{k_{m}J_{p}}{R_{m}J_{t}}\\\
\frac{k_{m}m_{p}lr}{R_{m}J_{t}}\end{array}\right],\;B_{1}={\rm
diag}\\{0.012,0.012,1,1\\},$ $\displaystyle C_{2}$
$\displaystyle=\left[\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\end{array}\right],\;D_{21}=\left[\begin{array}[]{cccc}0&0&0.0001&0\\\
0&0&0&0.0001\end{array}\right].$

where $J_{t}=J_{p}J_{r}-m_{p}^{2}l^{2}r^{2}$. The values of system parameters
of the rotary pendulum are presented in Table I. Since $w$ is considered as
linearization errors and deterministic unmodelled disturbances, the
disturbance matrices $B_{1}$ and $D_{21}$ above are chosen by trial and error.
Let

$\displaystyle\tilde{y}$
$\displaystyle=\left[\begin{array}[]{cc}1&0\end{array}\right]y,$ (24)

such that $E=\left[\begin{array}[]{cccc}1&0&0&0\end{array}\right]$, meaning
that the base arm angle $\theta_{1}$ needs to track a reference signal and the
vertical pendulum needs to be balanced in the upright position.

TABLE I: System parameters of the rotary pendulum. Parameter | Description | Value | Unit
---|---|---|---
$R_{m}$ | Terminal resistance | 8.4 | $\Omega$
$k_{m}$ | Motor back-emf constant | 0.042 |  V$\cdot$s/rad
$r$ | Length of the base arm | 0.085 | m
$2l$ | Length of the pendulum | 0.129 | m
$m_{p}$ | Mass of the pendulum | 0.024 | Kg
$g$ | Gravity constant | 9.81 | N/Kg
$b_{r}$ | Damping on the rotary arm | 0.0005 |  N$\cdot{\rm m}\cdot$ s/rad
$b_{p}$ | Damping on the pendulum | 0.0001 |  N$\cdot{\rm m}\cdot$ s/rad
$J_{r}$ | Inertia of the rotary arm | $2.3060\times 10^{-4}$ |  N$\cdot{\rm m}\cdot{\rm s}^{2}$/rad
$J_{p}$ | Inertia of the pendulum | $1.3313\times 10^{-4}$ |  N$\cdot{\rm m}\cdot{\rm s}^{2}$/rad

The weighting matrices for the performance variable $z$ in (6) and for the
tracking performance (15) are selected as
$C_{1}=\left[\begin{array}[]{cc}15&0\\\
\end{array}\right]^{\prime},D_{12}=\left[\begin{array}[]{cccc}0&\sqrt{2}\end{array}\right]^{\prime},R=C_{1}^{\prime}C_{1}$
and $Q=D_{12}^{\prime}D_{12}$. Then according to Step 2, we need to design the
LQT controller with $\tilde{Q}_{\alpha}$. The LQT controller in Proposition 1
is first designed for the linearized pendulum system, in which the closed-loop
poles are located at ($-16.55\pm j12.80,-21.20\pm j1.76$). The observer gain
$L$ is chosen such that the eigenvalues of $A+LC_{2}$ are ($-59.40\pm
j80.54,-61.04\pm j76.24$). The filter $\tilde{Q}$ is then designed to satisfy
the $H_{\infty}$ performance (13) [8], where $\gamma=0.21$. Following Step 3,
we shall tune $\alpha$ via the ES algorithm (21) in simulation under different
disturbance signals, shown in the next subsection.

###  V-B Simulation results

We consider that the disturbance signal $w$ is injected into the control
voltage $u$ and the measurement output $y$ is perturbed by white noise. In
particular, two different disturbance signals are considered: a square wave
with the amplitude of $2$ and frequency of 0.5Hz, denoted by $w_{1}$, and a
kind of vanishing signal $w_{2}(t)=5e^{-0.1t}$. Two different references are
considered: Case A tracks a square wave with amplitude of 20 degrees and
frequency of 0.05Hz while Case B tracks a sinusoidal signal
$r=\frac{\pi}{3}\sin\pi t$. In the simulation, the initial value of ES is
selected as $\alpha(0)=1$, indicating the local $H_{\infty}$ performance at
the first iteration.

For all simulations, the quadratic cost function $J(\alpha)$ in (15) is used
with $T=20s$ and the parameters $a$, $h$ and $\beta$ in the ES algorithm (21)
are set to

$a=0.8,h=0.1,\beta=0.015.$

The choice of the step size $\delta$ in the ES algorithm depends on a specific
disturbance signal and reference signal. The optimal $\alpha$ tuned via ES is
denoted as $\alpha^{*}$.

Now we tune $\alpha$ for disturbances $w_{1}$ and $w_{2}$. Table II summarizes
the $\alpha^{*}$ values in which minimum costs are reached for Case A and Case
B in the presence of $w_{1}$ and $w_{2}$. Figs. 6 and 7 show that ES minimizes
the cost function (15) with convergence to the parameter $\alpha$ that
produces a minimum for Case A and Case B in the presence of $w_{1}$
respectively. It can be seen that the controller with $\alpha^{*}$ obtained
from ES tuning yields a much better closed-loop performance in terms of
transient behaviours compared with the cases when $\alpha=0$ (LQT performance)
and $\alpha=1$ ($H_{\infty}$ performance). Figs. 6 and 7 show the
effectiveness of the proposed ES algorithm when dealing with deterministic
uncertainties and linearization errors. Similar results are observed when
dealing with $w_{2}$.

TABLE II: $\alpha^{*}$ for disturbances $w_{1}$ and $w_{2}$. | $w_{1}$ | $w_{2}$
---|---|---
Case A | $\alpha^{*}=1.20$ | $\alpha^{*}=1.21$
Case B | $\alpha^{*}=1.38$ | $\alpha^{*}=1.47$
Figure 6: Case A: ES $\alpha$ tuning for $w_{1}$ illustrated by (a) the
evolution of the cost function and (b) the parameter $\alpha$ during ES tuning
of the closed-loop system with $w_{1}$. The lower plots present (c) the output
signal $\theta_{1}$ and (d) the control input signal $u$. Figure 7: Case B: ES
$\alpha$ tuning for $w_{1}$ illustrated by (a) the evolution of the cost
function and (b) the parameter $\alpha$ during ES tuning of the closed-loop
system with $w_{1}$. The lower plots present (c) the output signal
$\theta_{1}$ and (d) the control input signal $u$.

##  VI Conclusion

In this paper, we develop a novel tracking control scheme for a classe of
nonlinear systems in the presence of disturbances based on the robust
controller with a Youla-type filter $\tilde{Q}$ and the data-driven ES
technique. A key point is that a gain factor $\alpha$ is introduced to the
filter $\tilde{Q}$, which allows an extra degree of freedom to be optimized
for tracking performance. It has been shown from the simulation results that
the proposed tracking controller with a tuned $\alpha$ via ES can achieve an
optimal tracking performance for nonlinear systems with disturbances.

## References

  * [1] M. Aliyu, _Nonlinear ${H_{\infty}}$-control, Hamiltonian systems and Hamilton-Jacobi equations_. Boca Raton, FL: CRC Press, 2011.
  * [2] H. Modares, F. L. Lewis, and Z.-P. Jiang, “${H_{\infty}}$ tracking control of completely unknown continuous-time systems via off-policy reinforcement learning,” _IEEE Transactions on Neural Networks and Learning systems_ , vol. 26, no. 10, pp. 2550–2562, 2015.
  * [3] H. Elmali and N. Olgac, “Robust output tracking control of nonlinear mimo systems via sliding mode technique,” _Automatica_ , vol. 28, no. 1, pp. 145–151, 1992.
  * [4] H. K. Khalil, _Nonlinear systems_ , 3rd ed. Upper Saddle River, NJ: Prentice-Hall, 2002.
  * [5] K. Zhou, J. Doyle, and K. Glover, _Robust and optimal control_. NJ: Prentice Hall, 1996.
  * [6] K. Zhou and Z. Ren, “A new controller architecture for high performance, robust, and fault-tolerant control,” _IEEE Transactions on Automatic Control_ , vol. 46, no. 10, pp. 1613–1618, 2001.
  * [7] B. D. Anderson and J. B. Moore, _Optimal control: Linear quadratic methods_. Englewood Cliffs, NJ: Prentice Hall, 1989.
  * [8] X. Chen, K. Zhou, and Y. Tan, “Revisit of LQG control–A new paradigm with recovered robustness,” in _Porceedings of the IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 5819–5825.
  * [9] T. He, X. Chen, and G. G. Zhu, “A dual-loop robust control scheme with performance separation: Theory and experimental validation,” _IEEE Transactions on Industrial Electronics_ , vol. 69, no. 12, pp. 13 483–13 493, 2022.
  * [10] D. J. Limebeer, B. D. Anderson, and B. Hendel, “A Nash game approach to mixed $H_{2}$/$H_{\infty}$ control,” _IEEE Transactions on Automatic Control_ , vol. 39, no. 1, pp. 69–82, 1994.
  * [11] J. Doyle, K. Zhou, K. Glover, and B. Bodenheimer, “Mixed $H_{2}$ and $H_{\infty}$ performance objectives II: Optimal control,” _IEEE Transactions on Automatic Control_ , vol. 39, no. 8, pp. 1575–1587, 1994.
  * [12] X. Chen and K. Zhou, “Multiobjective ${H_{\infty}}$/${H_{2}}$ control design,” _SIAM Journal on Control and Optimization_ , vol. 40, no. 2, pp. 628–660, 2001.
  * [13] K. B. Ariyur and M. Krstic, _Real-time optimization by extremum-seeking control_. NJ: John Wiley & Sons, 2003\.
  * [14] N. J. Killingsworth and M. Krstic, “PID tuning using extremum seeking: Online, model-free performance optimization,” _IEEE Control Systems Magazine_ , vol. 26, no. 1, pp. 70–79, 2006.
  * [15] Y. Tan, W. H. Moase, C. Manzie, D. Nešić, and I. M. Mareels, “Extremum seeking from 1922 to 2010,” in _Proceedings of the 29th Chinese Control Conference (CCC)_ , 2010, pp. 14–26.
  * [16] Q. Tan, P. Divekar, Y. Tan, X. Chen, and M. Zheng, “Model-guided extremum seeking for diesel engine fuel injection optimization,” _IEEE/ASME Transactions on Mechatronics_ , vol. 23, no. 2, pp. 936–946, 2018.
  * [17] B. Hu, K. Zhang, N. Li, M. Mesbahi, M. Fazel, and T. Başar, “Towards a theoretical foundation of policy optimization for learning control policies,” _arXiv preprint arXiv:2210.04810_ , 2022.
  * [18] J. Doyle, K. Glover, P. Khargonekar, and B. Francis, “State-space solutions to standard ${H}_{2}$ and ${H}_{\infty}$ control problems,” _IEEE Transactions on Automatic Control_ , vol. 34, no. 8, pp. 831–847, 1989.
  * [19] J.-Y. Choi, M. Krstic, K. B. Ariyur, and J. S. Lee, “Extremum seeking control for discrete-time systems,” _IEEE Transactions on Automatic Control_ , vol. 47, no. 2, pp. 318–323, 2002.
  * [20] J. Apkarian and M. Levis, _QUBE-Servo 2 Workbook (Student): Experiment for Matlab/Simulink Users, Quanser, Inc._ , 2021.