# Imaging the Dusty Substructures due to Terrestrial Planets in Planet-forming
Disks with ALMA and the Next Generation Very Large Array

Sarah K. Harter  Department of Physics and Astronomy, California State
University Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA  Luca
Ricci  Department of Physics and Astronomy, California State University
Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA  Shangjia Zhang
Department of Physics and Astronomy, University of Nevada Las Vegas, Las
Vagas, NV 89154, USA  Zhaohuan Zhu  Department of Physics and Astronomy,
University of Nevada Las Vegas, Las Vagas, NV 89154, USA Luca Ricci
<EMAIL_ADDRESS>
###### Abstract

We present simulations of the capabilities of the Atacama Large
Millimeter/submillimeter Array (ALMA) and of a Next Generation Very Large
Array (ngVLA) to detect and resolve substructures due to terrestrial planets
and Super-Earths in nearby planet-forming disks. We adopt the results of
global 2-D hydrodynamical planet-disk simulations that account for the
dynamics of gas and dust in a disk with an embedded planet. Our simulations
follow the combined evolution of gas and dust for several thousand planetary
orbits. We show that long integrations (several tens of hours) with the ngVLA
can detect and spatially resolve dust structures due to low-mass rocky planets
in the terrestrial planet formation regions of nearby disks (stellocentric
radii $r=1-3$ au), under the assumption that the disk viscosity in those
regions is low ($\alpha\leq 10^{-5}$). ALMA is instead unable to resolve these
structures in these disk regions. We also show that high-resolution ngVLA
observations separated by several days to few weeks would allow to detect the
proper motion of the azimuthally asymmetric structures expected in the disk
regions of terrestrial planet formation.

protoplanetary disks — circumstellar matter — planets and satellites:
formation

##  1 Introduction

The mutual interactions between planets and the parental disk play a prominent
role in setting the main physical properties of planetary systems. Planets are
made of solids and gas which come from the disk, and the mutual gravitational
interactions between a newly born planet and the material in the disk can
strongly affect the history of its location within the planetary system (for
reviews see Kley & Nelson, 2012; Baruteau et al., 2014).

Recent high angular resolution observations at both optical/near-infrared and
sub-millimeter/millimeter wavelengths have shown that substructures in disks
which resemble the predictions from planet-disk interactions theory are
ubiquitous in young nearby disks (e.g., Juhász et al., 2015; Dipierro et al.,
2018; Andrews et al., 2018b). The detailed modeling of these substructures,
especially rings and gaps, is often carried on via hydrodynamical simulations
to infer key properties of the putative young planets, such as mass and
orbital radius. These studies show that the rings and gaps seen by ALMA from
observations of the disk dust continuum likely trace a population of
relatively massive young planets, with masses higher than $10~{}M_{\oplus}$,
and relatively far from the host stars, at orbital radii greater than $10$ au
(Zhang et al., 2018; Long et al., 2018; Lodato et al., 2019).

Owing to different observational biases, the populations of these young
planets and of the detected mature exoplanets have a very scarce overlap in
terms of the main planetary and orbital properties, such as for example planet
mass and orbital radius (see Figure 21 in Zhang et al., 2018). This strongly
hinders comparisons between the properties of these young vs mature exoplanet
populations which are key to unveil the processes of formation and early
evolution of planetary systems (see discussion in Lodato et al., 2019). Future
improvements on the observational capabilities regarding both direct and
indirect methods of exoplanet detection will be necessary to inform theories
of planet formation and evolution.

Over the past several years, the astronomy community has shown strong interest
for a future large area mm/radio interferometric array optimized for imaging
at milliarcsecond (mas) resolution, which would be ideal for observations of
protoplanetary disks, among other sources in the sky. This Next Generation
Very Large Array (ngVLA, Murphy, 2018) would operate at frequencies of 1.2-116
GHz with approximately 10 $\times$ the collecting area of the current Very
Large Array (VLA) and Atacama Large Millimeter/submillimeter Array (ALMA)
interferometers, and with more than 10 $\times$ longer baselines that would
yield mas-resolution at the highest frequencies.

With the expected angular resolution and sensitivity at mm/cm wavelengths, the
ngVLA would be able to spatially resolve structures in the inner terrestrial
planet formation regions (1 au at the distance of the closest star forming
regions, i.e. 140 pc, corresponds to an angular scale of about 7 mas). In this
study we aim to quantify the capabilities of the ngVLA to detect and spatially
resolve disk structures due to the planets with masses lower than
$10~{}M_{\oplus}$ at orbital separations between 1 and 3 au from the host
star, which are more in line with the properties of the bulk of the population
of mature exoplanets detected so far.

Ricci et al. (2018) already tested the capabilities of the ngVLA to image disk
substructures due to planets in the dust continuum emission of nearby disks.
However, the simulations that were run for that study followed the evolution
of their disk-planet systems for less than 2000 planetary orbits. The models
presented here follow the disk evolution for a larger number of orbits, giving
more time for low-mass terrestrial planets to open gaps in the dust density of
low-viscosity disks, which can be detected by the ngVLA (see, e.g. Dong et
al., 2018).

Section 2 presents the methods used for this study, describing the disk
physical models used to simulate the disk-planet interactions, the procedure
used to obtain the model synthetic images at the ngVLA and ALMA wavelengths,
and the method used to simulate the results of future ngVLA and ALMA
observations. Section 3 outlines the results of this investigation, Section 4
presents a discussion of the main results, whereas Section 5 summarizes the
conclusions of this work.

##  2 Methods: from disk simulations to ngVLA and ALMA observations

We describe here the methods used to simulate observations with the ngVLA and
ALMA for the dust continuum emission of disks with embedded low-mass
terrestrial planets.

###  2.1 Disk model

The models used for this investigation were obtained via 2D hydrodynamical
planet–disk simulations using the modified version of the grid-based code
FARGO (Masset, 2000) called Dusty FARGO-ADSG (Baruteau & Masset, 2008a, b;
Baruteau & Zhu, 2016). This code accounts for the dynamics of gas and dust in
a disk with one or more embedded planets. In particular, it can calculate the
gravitational interaction between gas and dust with the planet, as well as the
aerodynamical gas-drag on the dust.

We refer to Zhang et al. (2018) for a detailed description of the disk models
and the method to extract synthetic images for the dust continuum at different
wavelengths for a given model. In this work, we initialized the disk structure
assuming azimuthal symmetry with a power-law radial dependence for the gas
surface density:

$\Sigma{{}_{g}}(r)=\Sigma_{g,0}(r/r{{}_{0}})^{-1},$ (1)

where $r$ is the stellocentric radius in the disk, and $r_{0}$ is an arbitrary
radius, that we assume to correspond to the orbital radius $r_{p}$ of the
planet in the simulations. The normalization factor for the gas surface
density $\Sigma_{g,0}$ is an input parameter for our models, and the values
considered for the initial conditions of the hydrodynamical simulations
presented here are $\Sigma_{g,0}=100,300,600,1200,2400$ g/cm2. Observations of
young disks conducted so far do not provide firm constraints on the gas
surface density at stellocentric radii $r<5-10$ au. The relatively large range
of $\Sigma_{g,0}$ values sampled in this work reflects the large spread
estimated for dust disk masses and disk radii in nearby star forming regions
(e.g., Andrews et al., 2013, 2018a). For reference, our largest value of
$\Sigma_{g,0}$ is a factor of about 1.4 greater than the gas surface density
at 1 au in the Minimum Mass Solar Nebula (MMSN) disk model presented in
Hayashi (1981), and a factor of 4 lower than in the MMSN model in Desch
(2007).

These models also assume the locally isothermal equation of state, where the
temperature at stellocentric radius $r$ is given by
$T(r)=T{{}_{0}}(r/r{{}_{0}})^{-1/2}$. The disk temperature $T$ is related to
the disk pressure scale height $h$ via $h/r=c_{s}/v_{\phi}$, where
$c_{s}=\sqrt{RT/\mu}$ is the local sound speed of the gas ($\mu=2.35$ is the
adopted value for the mean molecular weight), and $v_{\phi}$ is the local
rotation velocity of the gas, assumed to be keplerian. The disk aspect ratio
at the orbital radius $r_{p}=r_{0}$ of the planet is assumed to be
$h/r~{}(r_{p})=0.03$, and the models consider a Shakura-Sunyaev viscosity with
an $\alpha$ parameter of $10^{-5}$ constant across the disk (Shakura &
Sunyaev, 1973).

For the luminosity of the central star we considered two values of
$L_{\star}=1$ and $10~{}L_{\odot}$. In these models, the stellar luminosity
$~{}L_{\star}$ determines the temperature across the disk via

$T_{d}(r)=\left(\dfrac{\phi L_{\star}}{8\pi r^{2}\sigma}\right)^{1/4},$ (2)

where $\sigma$ is the Stefan-Boltzmann constant, and $\phi$ is a factor, taken
to be 0.02 from D’Alessio et al. (2001), which quantifies the effect of the
flaring of the disk on its heating from the stellar light.

We assume that the initial dust surface density of the disk is $1/100$ of the
gas surface density, and the dust is characterized by a grain size
distribution $n(s)\propto s^{-3.5}$, $s$ being the radius of the grain, with
maximum grain size of 1 cm. Since our dust component uses open boundary
conditions, the dust-to-gas mass ratio for the whole disk can change with
time. Dust and gas experience both gravitational forces and the dust feels
aerodynamic drag forces induced by gas. The Stokes number for particles, which
controls the level of aerodynamical coupling, is given by
$\rm{St}=t_{\rm{stop}}\Omega=\pi s\rho_{p}/2\Sigma_{g}$, where $t_{\rm{stop}}$
is the stopping time, $\Omega$ the angular velocity of the gas, and $\rho_{p}$
the solid density of dust particles. The combined temporal evolution of dust
and gas in the hydrodynamical code adopted in this work is described in more
detail in Zhang et al. (2018).

The numerical grid extends from $0.1r_{0}$ to $10r_{0}$ in the radial
direction and from 0 to $2\pi$ in the azimuthal direction, and is sampled by
750 and 1024 grid points in the radial and azimuthal directions, respectively.

In these models, the parameter $q$ represents the ratio between the planet
mass $M_{p}$ and the stellar mass $M_{\star}$ (see Zhang et al., 2018). We
have investigated two families of models, one with
$q=M_{p}/M_{\star}=1~{}M_{\oplus}/M_{\odot}$ and with
$10~{}M_{\oplus}/M_{\odot}$. The first family corresponds to having an Earth-
mass planet orbiting a Solar-mass star, or a Mars-mass planet around a M-type
star with mass of $0.1~{}M_{\odot}$. The second family corresponds to having a
$10~{}M_{\oplus}$-mass Super Earth around a Solar-mass star, or an Earth-mass
planet around a star with mass of $0.1~{}M_{\odot}$.

The simulations with $q=1~{}M_{\oplus}/M_{\odot}$ were run for up to 7000
planet orbits with the planet located at either $r_{p}=$ 2 or 3 au from the
host star. For the simulations with $q=10~{}M_{\oplus}/M_{\odot}$, the
simulations were run for up to 5000 planet orbits with the planet at $r_{p}=$
1 or 3 au. A summary of the values for the model parameters investigated in
this study is presented in Table 1.

The synthetic maps for the dust continuum emission of the disk models were
derived using the methods described in Zhang et al. (2018), and were produced
at wavelengths of 1.25, 3, 7 mm and 1 cm, and two representative cases are
shown in Figure 1.

$q$ [$M_{\oplus}/M_{\odot}$] | Planet orbital radius [au]
---|---
| $L_{\star}=10~{}L_{\odot}$ | $L_{\star}=1~{}L_{\odot}$
| $\Sigma_{g,0}^{\rm{init}}[\rm{g/cm}^{2}]=100$ | 300 | 600 | 1200 | 2400 | $\Sigma_{g,0}^{\rm{init}}=300~{}\rm{g/cm}^{2}$
1 | .. | 2, 3 | .. | .. | .. | 3
10 | 3 | 1, 3 | 3 | 3 | 3 | 1
Table 1: Values of the model parameters investigated in this study. Figure 1:
Model synthetic images for the dust continuum emission at 1.25, 3, 7 and 10
mm, respectively, for two of the systems considered in this study. Top row)
Disk model with a planet at 3 au from the host star and a planet-to-star mass
ratio of $1~{}M_{\oplus}/M_{\odot}$ at 7000 planet orbits as discussed in
Section 3.2. The planet location in each panel is marked with a green cross.
The integrated flux densities for this model are 42.0, 4.2, 0.32 and 0.099 mJy
at the wavelengths of 1.25, 3, 7 and 10 mm, respectively. Bottom row) Disk
model with a planet at 3 au from the host star and a planet-star mass ratio of
$10~{}M_{\oplus}/M_{\odot}$ at 5000 planet orbits as discussed in Section 3.1.
The integrated flux densities for this model are 18.7, 2.05, 0.170 and 0.055
mJy at the wavelengths of 1.25, 3, 7 and 10 mm, respectively.

###  2.2 Simulations of the ngVLA and ALMA observations

The synthetic maps obtained according to the methods outlined in Section 2.1
were converted into predictions for future observations with the ngVLA (at
$\lambda=3,7,10$ mm) and ALMA ($\lambda=1.25$ mm) using the CASA software
package111https://casa.nrao.edu/. We adopted the same procedure as in Ricci et
al. (2018), in which the ngVLA simulations are performed using the SIMOBSERVE
task to generate the visibility dataset in the $(u,v)$ Fourier space, and the
SIMNOISE task to add the noise by corrupting the visibilities. According to
the specifics of the reference design described in Selina et al. (2018), the
rms noise levels on the ngVLA maps presented in this paper correspond to on-
source integrations times of $\approx 70-80$ hours.

For the ngVLA simulations we considered the original ngVLA Rev B array
configuration distributed across the US Southwest and Mexico. This array
configuration includes 214 antennas of 18-meter diameter, to baselines up to
1000 km (Selina et al., 2018). In our simulations we did not include the 30
additional antennas distributed to baselines up to 9000 km, as the signal from
the sources considered here would be very faint for these very long baselines.
For the array configuration of the ALMA observations, we used the
alma.out28.cfg antenna position file available in the CASA package, which
incorporates the longest 16 km baselines in the ALMA array.

For the imaging of the visibilities we employed the CLEAN algorithm with
Briggs weighting, and adjusted the robust parameter to give a reasonable
synthesized beam and noise performance. In particular, the ALMA images were
computed with a Briggs weighting with robust parameter $R=-2$ (uniform
weighting), while for the ngVLA we chose $R=-1$. For the deconvolution we
employed a multiscale clean approach with scales of 0, 7 and 30 pixels to
better recover compact emission at both high brightness and large diffuse
structures in the disk models. The disk center was located at a declination of
$+35.0$ degrees for the ngVLA simulations and $-24.0$ degrees for the ALMA
simulations, and the assumed distance is 140 pc, similar to the distance of
nearby star forming regions such as Taurus or Ophiuchus.

##  3 Results

To facilitate the discussion of the results of the simulations with different
parameters for the disk, planet and star, we start by first presenting the
results of one specific model in Section 3.1. In the following sections, we
discuss separately the results when changing the planet-to-star mass ratio,
disk surface density, planet orbital radius, and stellar luminosity,
respectively.

###  3.1 Reference model

Figure 2: The azimuthally averaged radial profiles of the gas (top row) and
dust density (middle) and dust-to-gas mass ratio (bottom) for models with a
planet-to-star mass ratio $q=10~{}M_{\oplus}/M_{\odot}$ (left) and
$1~{}M_{\oplus}/M_{\odot}$ (right) with a planet located at 3 au, and an
initial gas surface density of 300 g/cm2 at the location of the planet.
Different radial profiles for different snapshots during the simulations are
shown in different colors as labeled on each panel. In each panel, the
abscissa of the green cross represents the stellocentric distance of the
planet, while its ordinate corresponds to the value on the y-axis at the
beginning of the simulations.  Figure 3: ALMA and ngVLA simulated observations
for the dust continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm
(ngVLA) for the disk model with a planet at 3 au from the host star and
planet-to-star mass ratio of $10~{}M_{\oplus}/M_{\odot}$, stellar luminosity
$L_{\star}=10~{}L_{\odot}$, and disk with initial surface gas density of 300
g/cm2 at the location of the planet. This is the reference model introduced in
Section 3.1. Each column represents the specific snapshot of the simulation in
terms of number of planet orbits, from left to right of 1000, 2000, 3000,
4000, and 5000 orbits, respectively. Each row represents simulations of either
ALMA or ngVLA observations at the wavelengths specified in the figure. The
resulting FWHM size of the synthesized beams are 18 mas (2.5 au) $\times$ 12
mas (1.7 au) at a wavelength of 1.25 mm, 1.7 mas (0.23 au) $\times$ 1.3 mas
(0.19 au) at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm, and
4.2 mas (0.59 au) $\times$ 3.7 mas (0.52 au) at 10 mm, respectively. The
resulting rms noise on the maps are 2.8 $\mu$Jy beam-1 at 1.25 mm, 86 nJy
beam-1 at 3 mm, 35 nJy beam-1 at 7 mm, and 23 nJy beam-1 at 10 mm. The disk
model has integrated fluxes of 28.3 mJy at 1.25 mm, 6.6 mJy at 3 mm, 0.53 mJy
at 7 mm, and 0.16 mJy at 10 mm, respectively, calculated at the 5000 orbits
snapshot. Figure 4: ALMA and ngVLA simulated observations for the dust
continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm (ngVLA) for the
disk models with a planet with planet-to-star mass ratio of
$1~{}M_{\oplus}/M_{\odot}$ at 3 au from the host star, stellar luminosity
$L_{\star}=10~{}L_{\odot}$, and initial gas surface density of 300 g/cm2 at
the location of the planet. This is the model introduced in Section 3.2. Each
column represents the specific snapshot of the simulation in terms of planet
orbits, from left to right, of 1000, 2000, 4000, 6000, and 7000 orbits. The
resulting synthesized beams have FWHM $=$ 18 mas (2.5 au) $\times$ 12 mas (1.7
au) at 1.25 mm, 1.7 mas (0.23 au) $\times$ 1.3 mas (0.19 au) at 3 mm, 3.4 mas
(0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm, and 4.2 mas (0.59 au) $\times$
3.7 mas (0.52 au) at 10 mm. The resulting rms noise on the maps are 2800 nJy
beam-1 at 1.25 mm, 86 nJy beam-1 at 3 mm, 35 nJy beam-1 at 7 mm, and 23 nJy
beam-1 at 10 mm. The disk model has integrated fluxes of 29.0 mJy at 1.25 mm,
4.70 mJy at 3 mm, 0.523 mJy at 7 mm, and 0.157 mJy at 10 mm, respectively,
calculated at the 7000 orbits snapshot.  Figure 5: Temporal evolution for the
$\Delta$ parameter, defined as the ratio between the average dust density in
the gap at the orbital radius of the planet and the average density of the
peaks of the two closest rings from the planet. Models are with
$q=10~{}M_{\oplus}/M_{\odot}$ (left panel) and $1~{}M_{\oplus}/M_{\odot}$
(right) with planet at an orbital radius of 3 au, and for several initial gas
surface densities at the stellocentric radius of the planet, as labeled in
each panel. Figure 6: ALMA and ngVLA simulated observations for the dust
continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm (ngVLA) for the
disk model with a planet with planet-to-star mass ratio of
$10~{}M_{\oplus}/M_{\odot}$ at 3 au from the host star, stellar luminosity
$L_{\star}=10~{}L_{\odot}$, and an initial gas surface gas density of 600
g/cm2 at the location of the planet. This model is discussed in Section 3.3.
Each column represents the specific orbit from left to right of 1000, 2000,
3000, 4000, and 5000 orbits. The resulting synthesized beam has a FWHM $=$ 18
mas (2.5 au) $\times$ 12 mas (1.7 au) for 1.25 mm, 1.7 mas (0.23 au) $\times$
1.3 mas (0.19 au) for 3 mm wavelength, 3.4 mas (0.47 au) $\times$ 2.7 mas
(0.38 au) for 7 mm wavelength, and 4.2 mas (0.59 au) $\times$ 3.7 mas (0.52
au) for 10 mm wavelength at the assumed distance of 140pc for our disk model.
The resulting rms noise on the ngVLA map is 2800 nJy beam-1 for 1.25 mm
wavelength, 86 nJy beam-1 at 3 mm wavelength, 35 nJy beam-1 at 7 mm
wavelength, and 21 nJy beam-1 at 10 mm wavelength. With integrated fluxes of
37.1 mJy at 1.25 mm, 9.29 mJy for 3 mm, 0.806 mJy at 7 mm, and 0.258 mJy at 10
mm. Figure 7: ALMA and ngVLA simulated observations for the dust continuum
emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm (ngVLA) for disk models with
a planet with planet-to-star mass ratio of $10~{}M_{\oplus}/M_{\odot}$ at 3 au
from the host star, stellar luminosity $L_{\star}=10~{}L_{\odot}$, and an
initial gas surface density of 100 g/cm2 at the location of the planet. Each
column represents the specific snapshot of the simulation in terms of number
of planet orbits, from left to right of 1000, 2000, 3000, 4000, and 5000
orbits, respectively. The resulting synthesized beams have FWHM $=$ 18 mas
(2.5 au) $\times$ 12 mas (1.7 au) at 1.25 mm, 1.7 mas (0.23 au) $\times$ 1.3
mas (0.19 au) at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm,
and 4.2 mas (0.59 au) $\times$ 3.7 mas (0.52 au) at 10 mm. The resulting rms
noise levels on the maps are 2800 nJy beam-1 at 1.25 mm, 86 nJy beam-1 at 3
mm, 35 nJy beam-1 at 7 mm, and 21 nJy beam-1 at 10 mm. The disk model has
integrated fluxes of 23.2 mJy at 1.25 mm, 4.07 mJy at 3 mm, 0.280 mJy at 7 mm,
and 0.080 mJy at 10 mm, respectively, calculated at the 5000 orbits snapshot.
Figure 8: ALMA and ngVLA simulated observations for the dust continuum
emission at 3 mm, 7 mm, and 10 mm (ngVLA) for disk models with a planet with
planet-to-star mass ratio of $10~{}M_{\oplus}/M_{\odot}$ at 3 au from the host
star, stellar luminosity $L_{\star}=10~{}L_{\odot}$, at 5000 orbits. Each
column shows the results for different initial gas surface densities at the
planet orbital radius, as labeled on the top of each column. The resulting
synthesized beams have FWHM $=$ 1.7 mas (0.23 au) $\times$ 1.3 mas (0.19 au)
at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm, and 4.2 mas
(0.59 au) $\times$ 3.7 mas (0.52 au) at 10 mm. The resulting rms noise levels
on the maps are 2800 nJy beam-1 at 1.25 mm, 86 nJy beam-1 at 3 mm, 35 nJy
beam-1 at 7 mm, and 21 nJy beam-1 at 10 mm. The disk models have integrated
fluxes at 3 mm of 4.07 mJy, 6.6 mJy, 9.29 mJy, 12.2 mJy, and 14.1 mJy, at 7 mm
of 0.280 mJy, 0.53 mJy, 0.806 mJy, 1.01 mJy, and 1.42 mJy, at 10 mm of 0.08
mJy, 0.16 mJy, 0.258 mJy, 0.373 mJy, and 0.521 mJy for initial surface
densities between 100 g/cm2 and 2400 g/cm2, respectively.  Figure 9: ALMA and
ngVLA simulated observations for the dust continuum emission at 1.25 mm
(ALMA), 3 mm, 7 mm, and 10 mm (ngVLA), from left to right, for a disk model
with planet-to-star mass ratio of $1~{}M_{\oplus}/M_{\odot}$, planet at 2 au
from the host star, $L_{\star}=10~{}L_{\odot}$, and an initial gas surface gas
density of 300 g/cm2 at the planet orbital radius. This disk model is
discussed in Section 3.4. These images were produced after 7000 planet orbits
from the beginning of the simulations. The resulting synthesized beams have
FWHM $=$ 18 mas (2.5 au) $\times$ 12 mas (1.7 au) at 1.25 mm, 1.7 mas (0.23
au) $\times$ 1.3 mas (0.19 au) at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7 mas
(0.38 au) at 7 mm, and 4.2 mas (0.59 au) $\times$ 3.7 mas (0.52 au) at 10 mm.
The resulting rms noise levels on the maps are 2800 nJy beam-1 at 1.25 mm, 86
nJy beam-1 at 3 mm, 35 nJy beam-1 at 7 mm, and 21 nJy beam-1 at 10 mm. The
model integrated fluxes are 26.5 mJy at 1.25 mm, 4.25 mJy at 3 mm, 0.330 mJy
at 7 mm, and 0.105 mJy at 10 mm. Figure 10: ALMA and ngVLA simulated
observations for the dust continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm,
and 10 mm (ngVLA) for disk models with a planet at 1 au from the host star,
planet-to-star mass ratio of $10~{}M_{\oplus}/M_{\odot}$, stellar luminosity
$L_{\star}=10~{}L_{\odot}$, and an initial gas surface density of 300 g/cm2 at
the planet location. Each column shows the temporal evolution with snapshots,
from left to right, of 1000, 2000, 3000, 4000, and 5000 planet orbits from the
beginning of the simulations, respectively. The resulting synthesized beams
have FWHM $=$ 18 mas (2.5 au) $\times$ 12 mas (1.7 au) at 1.25 mm, 1.7 mas
(0.23 au) $\times$ 1.3 mas (0.19 au) at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7
mas (0.38 au) at 7 mm, and 4.2 mas (0.59 au) $\times$ 3.7 mas (0.52 au) at 10
mm. The resulting rms noise levels are 2800 nJy beam-1 at 1.25 mm, 86 nJy
beam-1 at 3 mm, 35 nJy beam-1 at 7 mm, and 21 nJy beam-1 at 10 mm. The model
integrated fluxes are 10.8 mJy at 1.25 mm, 1.62 mJy at 3 mm, 0.120 mJy at 7
mm, and 0.048 mJy at 10 mm, respectively, calculated at the 5000 orbits
snapshot. Figure 11: ALMA and ngVLA simulated observations for the dust
continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm (ngVLA) for disk
models with a planet at 1 au from the star, with planet-to-star mass ratio of
$10~{}M_{\oplus}/M_{\odot}$, stellar luminosity $L_{\star}=1~{}L_{\odot}$, and
an initial gas surface density of 300 g/cm2 at the planet location. Each
column represents the specific snapshot of the simulation in terms of number
of planet orbits, from left to right, of 1000, 2000, 3000, 4000, and 5000
orbits, respectively. The resulting synthesized beams have FWHM $=$ 18 mas
(2.5 au) $\times$ 12 mas (1.7 au) at 1.25 mm, 1.7 mas (0.23 au) $\times$ 1.3
mas (0.19 au) at 3 mm, 3.4 mas (0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm,
and 4.2 mas (0.59 au) $\times$ 3.7 mas (0.52 au) at 10 mm. The resulting rms
noise levels on the maps are 2800 nJy beam-1 at 1.25 mm, 86 nJy beam-1 at 3
mm, 35 nJy beam-1 at 7 mm, and 21 nJy beam-1 at 10 mm. The model has
integrated fluxes of 5.86 mJy at 1.25 mm, 1.02 mJy at 3 mm, 0.082 mJy at 7 mm,
and 0.025 mJy at 10 mm. Figure 12: ALMA and ngVLA simulated observations for
the dust continuum emission at 1.25 mm (ALMA), 3 mm, 7 mm, and 10 mm (ngVLA)
for the disk model with a planet with planet-to-star mass ratio of
$1~{}M_{\oplus}/M_{\odot}$ at 3 au from the host star, stellar luminosity
$L_{\star}=1~{}L_{\odot}$, and an initial gas surface density of 300 g/cm2 at
the planet location. This model is presented in Section 3.5. Each panel is for
the snapshot at 7000 planet orbits from the beginning of the simulation. The
resulting synthesized beams have FWHM $=$ 18 mas (2.5 au) $\times$ 12 mas (1.7
au) at 1.25 mm, 1.7 mas (0.23 au) $\times$ 1.3 mas (0.19 au) at 3 mm, 3.4 mas
(0.47 au) $\times$ 2.7 mas (0.38 au) at 7 mm, and 4.2 mas (0.59 au) $\times$
3.7 mas (0.52 au) at 10 mm. The resulting rms noise levels on the maps are
2800 nJy beam-1 at 1.25 mm, 86 nJy beam-1 at 3 mm, 35 nJy beam-1 at 7 mm, and
21 nJy beam-1 at 10 mm. The disk model have integrated fluxes of 15.1 mJy at
1.25 mm, 4.10 mJy at 3 mm, 0.386 mJy at 7 mm, and 0.120 mJy at 10 mm.

We start discussing the results of our simulations for one specific model,
which, for simplicity, we will refer to as the reference model. This was
calculated with a planet at an orbital radius of 3 au and with a planet-to-
star mass ratio of $10~{}M_{\oplus}/M_{\odot}$, in a disk with an initial gas
surface density of 300 g/cm2 at the orbital radius of the planet. The star in
this reference model has a stellar luminosity $L_{\star}=10~{}L_{\odot}$.

The left column of Figure 2 shows the time evolution of the azimuthally-
averaged radial profiles of the gas and dust surface densities, and of the
local dust-to-gas mass ratio for this reference model. In particular, the top
panel shows how the effects of the disk-planet interaction are evident already
after 1000 planet orbits, and some evolution of the gas density is still
present at $4000-5000$ orbits. As already reported in several previous studies
for the disk-planet interaction in low-viscosity disks (e.g., Zhu et al.,
2014; Bae et al., 2016, 2017; Dong et al., 2017; Miranda & Rafikov, 2019),
multiple local maxima and minima are found in the radial profile of the gas
density outside of the location of the main gap centered at the planet orbital
radius, with a magnitude of these fluctuations strongly dependent on the
distance from the planet.

As a consequence of the aerodynamic coupling between gas and dust, and the
related radial drift of grains due to the drag by the gas rotating at non-
Keplerian speeds (e.g., Weidenschilling, 1977; Brauer et al., 2008), the
perturbations on the densities of dust have significantly higher amplitudes
than those for the gas. This results in multiple dusty rings which become more
and more prominent over time. Another consequence is that the local dust-to-
gas mass ratio oscillates between values well above 0.01 in the rings, down to
values well below 0.01 in the gaps. For reference, the Stokes parameter of
grains with sizes of 1 mm and 10 mm are of the order of $10^{-4}$ and
$10^{-3}$, respectively, at the beginning of the simulation for the reference
model.

Although the radial profiles in Figure 2 are useful to visualize and quantify
the effects of the planet-disk interactions on both the gas and dust
components of the disk, it is important to notice that the reference model
discussed in this section contains azimuthally asymmetric structures. These
are particularly evident in the images of the dust continuum emission at 3, 7
and 10 mm (bottom row in Figure 1). These structures are due to dust drifting
toward local pressure and density maxima produced by Rossby Wave instabilities
(RWIs) at the gap edges (e.g., Lovelace et al., 1999; Li et al., 2000, 2001;
Lyra & Umurhan, 2019). The fact that these structures are better visible at
longer wavelengths reflects the more efficient drift for grains of larger
sizes, which dominate the emission at longer wavelengths, and make the
azimuthal extent of the arc-like structure narrower for the dust emission at
lower frequencies. Another important factor is the lower optical depth of the
dust emission at longer wavelengths which makes the surface brightness more
sensitive to the spatial variation of dust density, hence increasing the
contrast on the map. For example, in the reference model after 5000 planet
orbits, the optical depth of the dust emission near the peaks of the rings is
$\approx 20$ at 1 mm and $\approx 1$ at 1 cm. For the model discussed in the
next section with $q=1~{}M_{\oplus}/M_{\odot}$ the optical depth is $\approx
7$ at 1 mm and $\approx 0.4$ at 1 cm.

These structures are visible in Figure 3, which showcases the results of the
ALMA and ngVLA observations of the reference model for different snapshots of
the simulations. Whereas the gap structure is just barely visible in the ALMA
map at 1.25 mm, the ngVLA observations are able to detect and spatially
resolve several radial and azimuthal structures in the disk.

The comparison between maps for different snapshots also highlight a clear
evolutionary trend, i.e. nearly all the dusty structures tend to become more
and more pronounced with time, as a consequence of the radial drift of grains
discussed above.

For the first few thousand planet orbits, compact emission from dust in the
Lagrangian points L4 and L5 of the star-planet system is also detected (and
with best visibility at 3 mm), but already after 4000 orbits the emission gets
below the sensitivity limits of the simulated observations. These features are
also visible in the radial profile of the dust surface density shown in Figure
2. A depletion timescale of few thousand orbits for dust in these Lagrangian
points can have important consequences for the Trojan dust and larger Trojan
bodies expected for terrestrial planets and Super Earths (see Montesinos et
al., 2020, for a discussion of Trojan bodies of giant planets, and possible
implications for the Trojan solids found in the Solar System).

###  3.2 Models with lower planet-to-star mass ratio

After the presentation of the reference model, the first parameter we vary is
the planet-to-star mass ratio. The right column of Figure 2 shows the radial
profiles of the gas and dust density across the disk for a disk model with the
same parameters as the reference model, except for a planet-to-star mass ratio
10 times lower, i.e. $q=1~{}M_{\oplus}/M_{\odot}$, while keeping the planet
orbital radius at 3 au from the star. Because of the much weaker gravity of
the planet, the gravitational torques on the disk are weaker as well, and the
rings/gaps structures are much less pronounced than in the reference model,
both in terms of the contrast between rings and gaps and radial width of the
gaps.

At the same time, as for the reference model, the radial drift experienced by
grains significantly amplifies the structures in the dust component. For
example, after 7000 planet orbits, the dust surface density at the location of
the planet is about $1~{}\rm{g/cm}^{2}$, whereas towards the peaks in the
inner and outer rings it reaches $\approx 8-10~{}\rm{g/cm}^{2}$. In the gap,
the dust-to-gas mass ratio goes down to 0.003, while in the rings it reaches
values of $0.02-0.03$, hence with oscillations of about $2-3\times$ below and
above the initial ratio of 0.01.

Another important difference with the reference model is that the less steep
radial gradient in the gas density does not produce RWIs at the gap edges,
and, as a consequence, the disk structures are much more azimuthally
symmetric. This can be seen in Figure 4, which showcases the results of the
ALMA and ngVLA observations for the disk model with
$q=1~{}M_{\oplus}/M_{\odot}$. For this model, the gap opened by the planet is
detected by the ngVLA, with decreasing signal-to-noise ratios at longer
wavelengths, similarly to the reference model, whereas ALMA fails to detect
the gap as a consequence of the much larger synthesized beam. In the case of
the ngVLA map at 3 mm and at 7000 orbits, the peak signal-to-noise ratio on
the inward and outward edges of the gap itself are 12.1 and 11.2,
respectively.

We also notice that the ring/gap structure is not visible within $1000-2000$
planet orbits from the beginning of the simulations. Only after 2000 planet
orbits the gap starts to be marginally detected with the ngVLA at 3 mm, and
the detection becomes more and more significant after $3000-4000$ orbits.

###  3.3 Models with different gas surface densities

We present here the results of the models in which we varied the initial gas
surface density of the disk. With a fixed initial dust-to-gas mass ratio of
100, varying the gas density automatically modifies also the dust density.
Moreover, varying the gas density also changes the aerodynamic coupling
between the gas and dust grains (see above the definition of the Stokes
parameter), with an influence on the timescales of radial drift of solids
across the disk. The reference model presented in Section 3.1 considered an
initial gas surface density of 300 g/cm2 at the planet location ($r_{p}=$ 3
au). Here we discuss the results of simulations with values of 100, 600, 1200,
and 2400 g/cm2 for this parameter.

Figure 5 displays the $\Delta$ parameter, which quantifies the amount of dust
within the gap relative to the amount of dust in the rings, as a function of
time. The plots show the results of simulations with different values of the
initial gas surface density, and for values of the planet-to-star mass ratio
of $1~{}M_{\oplus}/M_{\odot}$ (left plot) and $10~{}M_{\oplus}/M_{\odot}$
(right). More precisely, the $\Delta$ parameter is defined as the ratio
between the average dust density in the gap at the orbital radius of the
planet and the average density of the peaks of the two closest rings from the
planet, one inside and one outside of the planet orbit. With this definition,
$\Delta$ depends both on the dust density in the gap and in the inner and
outer rings, but, as shown in Fig. 2, the variation of the density of dust
corotating with the planet is much larger than the density in the rings
(especially for the case with $10~{}M_{\oplus}/M_{\odot}$).

The comparison between the two panels in Figure 5 highlights the very
different timescales in the depletion of the gap in systems with different
planet-to-star mass ratios. In the simulations with
$q=10~{}M_{\oplus}/M_{\odot}$, the $\Delta$ parameter reaches values $\sim
0.2-0.7$ within the first $1000-3000$ planet orbits of evolution, and then
strongly decreases by factors of $\sim 10$ and $\sim 100$ at 4000 and 5000
orbits, respectively. Conversely, the temporal evolution is much more gradual
in the case of the simulations with $q=1~{}M_{\oplus}/M_{\odot}$. Here the
weaker gravity of the planet takes much longer to deplete the gap, as even
after 7000 planet $\Delta\sim 0.1-0.2$, and the decrease with time is roughly
linear. In terms of the dependence on the gas surface density of the disk,
although the qualitative trend does not depend on density, one can notice
faster evolution for the disk models with lower gas densities. This is in line
with the expectation of a faster radial drift of grains with larger (but lower
than 1) values of the Stokes parameter, since the Stokes parameter is
inversely proportional to the local density of gas in the disk.

The simulated observations for the disk models with initial gas surface
densities at 3 au of 600 and 100 g/cm2 are shown in Fig. 6 and 7,
respectively. The main difference between these results, including also the
results for the reference model with gas density of 300 g/cm2 shown in Figure
3, is that the models with lower gas densities have generally more pronounced
substructures: gaps are radially wider and with higher depletion, and the
azimuthally asymmetric structures related to local dust trapping are also more
concentrated along the azimuthal direction, including also in the Lagrangian
points L4 and L5 (Lyra et al., 2009). Similarly to the discussion above for
the depletion in the planet-induced gaps, this is an effect of the faster
drift of grains with higher values of the Stokes parameter. Despite the fact
that dust substructures tend to be more prominent in disk models with lower
gas densities, the correspondingly lower densities in dust make some of these
structures detectable only at marginal signal-to-noise ratios, especially at
10 mm (Fig. 7). At the same time, the higher optical depths for the models
with the highest densities considered in this work limit the visibility of the
azimuthal structures even at 10 mm.

###  3.4 Models with different planet orbital radii

Another model parameter that was varied in our simulations is the orbital
radius of the planet. This parameter has an effect on the substructures
generated from the disk-planet interactions. For example, the radial width of
the main gap increases linearly with the orbital radius of the planet (e.g.,
Ricci et al., 2018).

Whereas the reference model presented in Section 3.1 considered a planet at 3
au from the star, we also ran simulations for a disk model with
$q=1~{}M_{\oplus}/M_{\odot}$ with the planet at 2 au, and for a disk model
with $q=10~{}M_{\oplus}/M_{\odot}$ with the planet at 1 au.

Figure 9 displays the results of the ALMA and ngVLA simulations for models
with $q=1~{}M_{\oplus}/M_{\odot}$ and planet at 2 au. The maps show that the
gap is only detectable at wavelengths of 3 and 7 mm, barely detectable at 10
mm, and not detected with ALMA. The same gap is only marginally detected even
at 3 mm in the snapshots with planet orbits lower than 7000. We also ran
simulations for a planet with the same $q-$ratio at 1 au from the star, but
the gap was not detectable with neither ngVLA nor ALMA, and for that reason
those results are not shown here.

For the model with planet at 1 au and $q=10~{}M_{\oplus}/M_{\odot}$, Figure 10
showcases the temporal evolution of the dust emission at wavelengths of 1.25
mm, 3 mm, 7 mm, and 10 mm, respectively. As for the previous case discussed in
this section, the ALMA observations lack the angular resolution to detect the
gap opened by the planet. The situation is different for the ngVLA, which can
instead detect both the main gap structure as well as azimuthally asimmetric
structures, especially at 3 mm. The gap structure and the arc-like structure
is detected, although at lower signal-to-noise ratios, also at 7 and 10 mm
with the ngVLA, with the visibility of these structures significantly
increasing after few thousand planet orbits of evolution.

###  3.5 Models with lower stellar luminosity

In all the simulations presented so far, the bolometric luminosity of the star
was assumed to be $L_{\star}=10~{}L_{\odot}$. For several pre-Main Sequence
evolutionary models, a young star with the mass of the Sun reaches this
luminosity within its first Myr of age (e.g., Siess et al., 2000). In this
section we present the results of simulations with a lower stellar luminosity
of $L_{\star}=1~{}L_{\odot}$. In this work, this change affects only the dust
temperature which is used for the computation of the model images for the dust
continuum emission, but not the gas temperature, and pressure, used in the
hydrodynamical simulations. We also note that, at a fixed stellar mass, the
dependence of the vertical scale height with stellar luminosity is rather
weak. For example, a disk around a Solar-mass star with
$L_{\star}=1~{}L_{\odot}$ has an aspect ratio $h/r=0.03$ at 3 au assuming
$\phi=0.02$ (Eq. 2). If $L_{\star}=10~{}L_{\odot}$, then $h/r=0.04$. Hence,
the model images with lower stellar luminosity are characterized by a lower
surface brightness as a consequence of the lower dust temperature.

Figure 11 displays the results of the model with $q=10~{}M_{\oplus}/M_{\odot}$
and with planet at 1 au from a star with $L_{\star}=1~{}L_{\odot}$. These
simulations are similar to the case with $L_{\star}=10~{}L_{\odot}$ shown in
Fig. 10, but have lower surface brightnesses by a factor of $\approx 2$. Even
with these lower fluxes, the gap is still detectable with the ngVLA, with the
usual trends of higher signal-to-noise ratios at 3 mm and towards the later
snapshots of the simulations.

Figure 12 shows the results of the simulations after 7000 planet orbits for a
planet at 3 au from the star with $q=1~{}M_{\oplus}/M_{\odot}$, and again with
stellar luminosity $L_{\star}=1~{}L_{\odot}$. The planet-induced gap is
noticeable with the ngVLA at 3 mm and at 7 mm, even though with marginal
signal-to-noise ratio.

##  4 Discussion

Figure 13: Map for the spectral index of the dust emission between 3 and 7 mm
for the reference model presented in Section 3.1. This was derived by
combining the simulated ngVLA images at 3 and 7 mm after 5000 planet orbits,
shown in Figure 3.

In Section 3 we have shown that the ngVLA can resolve several substructures
due to the disk-planet interaction in low viscosity disk systems in nearby
star forming regions, with $q=1-10~{}M_{\oplus}/M_{\odot}$ and with the planet
at orbital radii down to $1-2$ au from the central star. These substructures
include multiple rings/gaps, arc-like structures and dust concentrated in
Lagrangian points of the planet-star system. Although the signal-to-noise
ratios decrease at longer wavelengths as a consequence of the coarser angular
resolution and weaker dust emission at lower frequencies, our results show
that for several models, especially those with $q=10~{}M_{\oplus}/M_{\odot}$,
the combined observations at multiple wavelengths, between 3 mm and 1 cm in
this study, can highlight variations in the morphology of some of these
structures across wavelength. This is a result of the spatial segregation of
grains with different sizes, as expected by the mechanism of radial drift of
solids in a gaseous disk (e.g., Laibe et al., 2012; Laibe, 2014; Laibe et al.,
2014; Birnstiel et al., 2016).

The variation across wavelength of the disk emission can be diagnosed also
through the map of the spectral index $\alpha_{\rm{mm}}$ of the dust emission
at mm wavelengths ($F_{\nu}\propto\lambda^{\alpha_{\rm{mm}}}$). Because of the
expected anti-correlation between the absolute value of $\alpha_{\rm{mm}}$ and
grain size, at least for grain sizes larger than $\sim 0.1-1$ mm (e.g., Ricci
et al., 2010), one would expect to find the lowest absolute values on the
spectral index map in correspondence to the dust substructures formed via the
trapping of the larger grains, as discussed above (van der Marel et al., 2015;
Barge et al., 2017). This is confirmed in Figure 13, which shows absolute
values of $\alpha_{\rm{mm}}$ as low as $\approx 1.8-2.2$ at the locations of
the dust substructures found in the reference model presented in Section 3.1.

It is important to emphasize that for several of the models presented here,
especially the ones with the lower $q$-value of $1~{}M_{\oplus}/M_{\odot}$,
these substructures start to become detectable in the dust continuum maps only
after several thousand planet orbits. This demonstrates the importance of
simulations that span at least several thousand orbital timescales of
evolution, especially when investigating the disk-planet interaction with low-
mass planets like Earth-like or Super Earth planets. This also explains the
apparent discrepancy between the results presented here and those from
previous similar investigations, such as for example in Ricci et al. (2018),
in which hydrodynamical simulations were run only for less than 2000 planet
orbits. This result is particularly interesting also considering that Super-
Earth and Earth-mass planets are more commonly found around Main Sequence
stars than higher mass giant planets (e.g., Howard et al., 2010; Petigura et
al., 2013).

The simulations of the observations presented in this work assumed a distance
of 140 pc for the disks. This is similar to the average distance of young
stars in nearby regions such as Taurus (e.g., Kenyon et al., 2008) or
Ophiuchus (e.g., Wilking et al., 2008). Andrews et al. (2013) performed a
survey in Taurus of 227 young stars surrounded by disks using the
Submillimeter Array (SMA) at a wavelength of 1.3 mm, which is comparable to
the wavelength of 1.25 mm for our ALMA simulations. 80 disks of those
surveyed, have a flux density at 1.3 mm above 10 mJy, which is close to the
flux density value of the faintest models in our investigation. An ALMA survey
of 147 disks was presented by Cieza et al. (2019) in the Ophiuchus star
forming region as part of the Ophiuchus Disc Survey Employing ALMA (ODISEA)
program. They found 60 disks with integrated fluxes greater than 10 mJy.

Another potentially interesting investigation that the ngVLA observations
would allow is the study of the proper motion of the azimuthally asymmetric
structures, which are nearly always detected in our simulations with
$q=10~{}M_{\oplus}/M_{\odot}$. For example, given the Keplerian orbital
velocity at 1 au from a Solar-mass star, and a spatial resolution of 0.2 au
from the ngVLA at 3 mm for a disk at 140 pc, the proper motion of a source
detected at a signal-to-noise ratio of $\approx 5$, would be detectable at
3$\sigma$ with two epochs separated by just $3-4$ days. For a source orbiting
at 3 au from a Solar-mass star the proper motion would be detectable at
3$\sigma$ with epochs separated by about 6 days.

We also note that our simulations lack the numerical resolution to investigate
the possible presence and emission of circumplanetary disks accreting onto the
planets (e.g., Isella et al., 2019). The potential for detection of these
systems by ALMA and the ngVLA has been discussed in Zhu et al. (2018).

Finally, although in this work we have investigated the dependence of the
results on a variety of different model parameters, there are other parameters
that we kept fixed for our investigation. For example, we fixed the $\alpha$
viscosity parameter to a value of $10^{-5}$. As shown in previous works (e.g.,
Fung et al., 2014; Kanagawa et al., 2016; Dong et al., 2017; Ricci et al.,
2018), gas with higher viscosity by factors $>10-100$ than considered in this
paper would make the density perturbations created by the interaction with the
planet to be more smoothed out, which would have an impact on the observed
structures in the dust emission. For the specific case of the disk models
considered here, Zhang et al. (2018) showed that in models with higher values
for the gas viscosity $\alpha=10^{-3}$ and $10^{-4}$, gaps are opened in the
case of $q\approx 10~{}M_{\oplus}/M_{\odot}$ (see their Fig. 6 and 7), but not
for $q\approx 1~{}M_{\oplus}/M_{\odot}$.

The simulations use a locally-isothermal equation of state. Recent works by
Miranda & Rafikov (2019, 2020) and Zhang & Zhu (2020) find that cooling can
affect the shape of the gaps. In the inner disk regions considered in this
study, the cooling timescale can be longer than the orbital timescale. Under
these conditions, considering an adiabatic equation of state can have
significant effects on the results of this work, and this will be the subject
of a follow-up work in the future.

Another parameter that we kept fixed in our simulations is the vertical aspect
ratio of the disk: $h/r~{}(r_{p})=0.03$. The vertical scale height sets the
strength of the pressure forces in the gas component of the disk, which,
similarly to viscosity, act against the formation of density perturbations in
the disk. We also did not investigate the dependence of our results on the
grain size distribution of dust in the disk. In our simulations we assumed an
initial grain size distribution $n(s)\propto s^{-3.5}$ with a maximum grain
size of 1 cm (see Zhang et al., 2018). We postpone the study on the impact of
these assumptions to a future work.

##  5 Conclusions

In this work we have investigated and tested the imaging capabilities of the
ngVLA, and the comparison with ALMA, for resolving gaps and rings and other
substructures in the dust continuum of disks with low viscosity due to low-
mass terrestrial planets. We produced models with
$q=1~{}M_{\oplus}/M_{\odot}$, corresponding for example to an Earth-mass
planet around a Solar-mass star or a Mars-mass planet around a
$0.1~{}M_{\odot}$ M-type star, and with $q=10~{}M_{\oplus}/M_{\odot}$,
corresponding to a Super-Earth around a Solar-mass star or an Earth-mass
planet around a 0.1 $M_{\odot}$ M-type star 222At the same time, it is
important to remind the reader that our simulations were run with a fixed
radial profile of the vertical scale height which is more consistent with the
observed properties of disks around young stars with masses $\approx
0.5-1.0~{}M_{\odot}$..

We showed that after several thousand orbits of evolution in these systems,
terrestrial planets with orbital radii of $1-3$ au can leave imprints in the
dust continuum emission that the ngVLA would be able to detect. Deep
observations, especially at 3 mm, would allow us to well characterize the
radial structure of the main rings/gaps which are found in several disk models
with $q=1$ and $10~{}M_{\oplus}/M_{\odot}$, and also detect emission from
azimuthal structures close to the planet for models with
$q=10~{}M_{\oplus}/M_{\odot}$ which are predicted by the physics of the
planet-disk interaction in disks with low viscosity.

Furthermore, multi-wavelength observations within the range of frequencies
covered by the ngVLA would allow us to map the emission dominated by grains
with different sizes within regions of local dust concentration. This would
provide further observational constraints to the spatial segregation of grains
with different sizes, to the measurement of the dust optical depth at these
wavelengths, and finally to further constrain the physical process of planet-
disk interaction and of grain growth in high-density regions which could be
suitable for the formation of planetesimals.

Given the very high angular resolution and astrometric precision expected for
the ngVLA, together with the short orbital timescales at few au from young
stars, the proper motion due to the orbital motion expected for these
structure should be detectable already with observations separated by few
days. This would be important also to confirm the physical association of
compact structures to the disk, rather than to background (or foreground)
sources which may be expected at the very high sensitivities probed by the
ngVLA.

Our results show that a ngVLA with the design considered in recent studies has
the potential to transform the fields of planet formation and protoplanetary
disks, allowing for imaging of the disk physical structures due to low-mass
terrestrial planets at sub-au resolution in nearby star forming regions.

Although this work presents results over a relatively broad range of the
parameter space of our models (planet-to-star mass ratio, planet orbital
radius, disk density, stellar luminosity), there are other assumptions that
can have an important impact on the results discussed here (e.g., gas
viscosity, grain size distribution, thermodynamics). Further investigations
over a broader region of the model parameter space will be provided in a
future work.

We thank the anonymous referee for his/her helpful comments. This work was
supported in part by the ngVLA Community Studies program, coordinated by the
National Radio Astronomy Observatory, which is a facility of the National
Science Foundation operated under cooperative agreement by Associated
Universities, Inc.

## References

  * Andrews et al. (2013) Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129, doi: 10.1088/0004-637X/771/2/129
  * Andrews et al. (2018a) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018a, ApJ, 865, 157, doi: 10.3847/1538-4357/aadd9f
  * Andrews et al. (2018b) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018b, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Bae et al. (2016) Bae, J., Nelson, R. P., & Hartmann, L. 2016, ApJ, 833, 126, doi: 10.3847/1538-4357/833/2/126
  * Bae et al. (2017) Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201, doi: 10.3847/1538-4357/aa9705
  * Barge et al. (2017) Barge, P., Ricci, L., Carilli, C. L., & Previn-Ratnasingam, R. 2017, A&A, 605, A122, doi: 10.1051/0004-6361/201629918
  * Baruteau & Masset (2008a) Baruteau, C., & Masset, F. 2008a, ApJ, 672, 1054, doi: 10.1086/523667
  * Baruteau & Masset (2008b) —. 2008b, ApJ, 678, 483, doi: 10.1086/529487
  * Baruteau & Zhu (2016) Baruteau, C., & Zhu, Z. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3927, doi: 10.1093/mnras/stv2527
  * Baruteau et al. (2014) Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 667
  * Birnstiel et al. (2016) Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41, doi: 10.1007/s11214-016-0256-1
  * Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859, doi: 10.1051/0004-6361:20077759
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698, doi: 10.1093/mnras/sty2653
  * D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, The Astrophysical Journal, 553, 321, doi: 10.1086/320655
  * Desch (2007) Desch, S. J. 2007, ApJ, 671, 878, doi: 10.1086/522825
  * Dipierro et al. (2018) Dipierro, G., Ricci, L., Pérez, L., et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, 5296, doi: 10.1093/mnras/sty181
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127, doi: 10.3847/1538-4357/aa72f2
  * Dong et al. (2018) —. 2018, ApJ, 866, 110, doi: 10.3847/1538-4357/aadadd
  * Fung et al. (2014) Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88, doi: 10.1088/0004-637X/782/2/88
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Howard et al. (2010) Howard, A. W., Marcy, G. W., Johnson, J. A., et al. 2010, Science, 330, 653, doi: 10.1126/science.1194854
  * Isella et al. (2019) Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25, doi: 10.3847/2041-8213/ab2a12
  * Juhász et al. (2015) Juhász, A., Benisty, M., Pohl, A., et al. 2015, Monthly Notices of the Royal Astronomical Society, 451, 1147, doi: 10.1093
  * Kanagawa et al. (2016) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, Publications of the Astronomical Society of Japan, 68, doi: 10.1093/pasj/psw037
  * Kenyon et al. (2008) Kenyon, S. J., Gómez, M., & Whitney, B. A. 2008, Low Mass Star Formation in the Taurus-Auriga Clouds, ed. B. Reipurth, Vol. 4, 405
  * Kley & Nelson (2012) Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211, doi: 10.1146/annurev-astro-081811-125523
  * Laibe (2014) Laibe, G. 2014, MNRAS, 437, 3037, doi: 10.1093/mnras/stt1928
  * Laibe et al. (2012) Laibe, G., Gonzalez, J. F., & Maddison, S. T. 2012, A&A, 537, A61, doi: 10.1051/0004-6361/201015349
  * Laibe et al. (2014) Laibe, G., Gonzalez, J.-F., & Maddison, S. T. 2014, MNRAS, 437, 3025, doi: 10.1093/mnras/stt1927
  * Li et al. (2001) Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874, doi: 10.1086/320241
  * Li et al. (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023, doi: 10.1086/308693
  * Lodato et al. (2019) Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, Monthly Notices of the Royal Astronomical Society, 486, 453, doi: 10.1093
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Lyra et al. (2009) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125, doi: 10.1051/0004-6361:200810797
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, Publications of the Astronomical Society of the Pacific, 131, 072001, doi: 10.1088/1538-3873/aaf5ff
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165, doi: 10.1051/aas:2000116
  * Miranda & Rafikov (2019) Miranda, R., & Rafikov, R. R. 2019, ApJ, 878, L9, doi: 10.3847/2041-8213/ab22a7
  * Miranda & Rafikov (2020) —. 2020, ApJ, 892, 65, doi: 10.3847/1538-4357/ab791a
  * Montesinos et al. (2020) Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, arXiv e-prints, arXiv:2009.10768.  https://arxiv.org/abs/2009.10768
  * Murphy (2018) Murphy, E. 2018, Science with a Next Generation Very Large Array, Vol. 517
  * Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Sciences, 110, 19273–19278, doi: 10.1073/pnas.1319909110
  * Ricci et al. (2018) Ricci, L., Liu, S.-F., Isella, A., & Li, H. 2018, ApJ, 853, 110, doi: 10.3847/1538-4357/aaa546
  * Ricci et al. (2010) Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15, doi: 10.1051/0004-6361/200913403
  * Selina et al. (2018) Selina, R. J., Murphy, E. J., McKinnon, M., et al. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 15
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593.  https://arxiv.org/abs/astro-ph/0003477
  * van der Marel et al. (2015) van der Marel, N., Pinilla, P., Tobin, J., et al. 2015, ApJ, 810, L7, doi: 10.1088/2041-8205/810/1/L7
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57, doi: 10.1093/mnras/180.1.57
  * Wilking et al. (2008) Wilking, B. A., Gagné, M., & Allen, L. E. 2008, Star Formation in the $\rho$ Ophiuchi Molecular Cloud, ed. B. Reipurth, Vol. 5, 351
  * Zhang & Zhu (2020) Zhang, S., & Zhu, Z. 2020, MNRAS, 493, 2287, doi: 10.1093/mnras/staa404
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744
  * Zhu et al. (2018) Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850, doi: 10.1093/mnras/sty1503
  * Zhu et al. (2014) Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122, doi: 10.1088/0004-637X/785/2/122


# Relevance of long-range screening in Mott transition examined via a hydrogen
lattice

Zi-Jian Lang (gbsn郎子健)  Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Sudeshna Sen
Department of Physics, Indian Institute of Technology (ISM) Dhanbad, India,
826004  Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai
Jiao Tong University, Shanghai 200240, China  Pak Ki Henry Tsang  Department
of Physics and National High Magnetic Field Laboratory, Florida State
University, Tallahassee, Florida 32306, USA  Kristjan Haule  Center for
Materials Theory, Department of Physics & Astronomy, Rutgers University,
Piscataway, NJ, O8854, USA  Vladimir Dobrosavljević<EMAIL_ADDRESS>Department of Physics and National High Magnetic Field Laboratory, Florida
State University, Tallahassee, Florida 32306, USA  Wei Ku (bsmi顧威)
<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China  Ministry of Education
Key Laboratory of Artificial Structures and Quantum Control, Shanghai 200240,
China

###### Abstract

The Mott transition, a metal-insulator transition due to strong electronic
interaction, is observed in many materials without an accompanying change of
system symmetry. An important open question in Mott’s proposal is the role of
long-range screening, whose drastic change across the quantum phase transition
may self-consistently make the transition more abrupt, toward a first-order
one. Here we investigate this effect in a model system of hydrogen atoms in a
cubic lattice, using charge self-consistent dynamical mean-field theory that
incorporates approximately the long-range interaction within the density
functional treatment. We found that the system is well within the charge-
transfer regime and that the charge-transfer gap intimately related to the
Mott transition closes smoothly instead. This indicates that within this
approximation, the long-range screening does not play an essential role in
this prototypical example. We explain this finding by demonstrating
analytically that the obtained insulating phase in this typical case is driven
by short-range physics of large energy scale, and the transition is controlled
by the charge-transfer gap closing, both insensitive to long-range screening.

Mott transition, as an interaction-driven metal-insulator transition (MIT)
without assistance of symmetry change, continues to attract a great attention
in the past decades [1, 2, 3]. It describes the phase transition between a
metal and the Mott insulator widely observed in strongly correlated materials,
such as transition-metal oxides [4, 5, 6, 7, 8]. To date, Mott insulators have
become one of the most important platforms for the exploration of strongly
correlated physics such as unconventional superconductivity [4, 9] and
magnetism [10].

The unique characteristic of Mott insulators lies in their unexpected
insulating behavior. Traditional insulators have their Bloch orbitals
completely filled by electrons such that a finite energy scale is needed for
the electrons to be able to propagate without violating the Pauli principle.
Mott insulators, on the other hand, do not satisfy this condition and
therefore would appear metallic in standard band theories, indicating the need
to include additional physics beyond the Pauli principle. To understand such
insulating behavior, Mott [11, 1] firstly pointed out the importance of
electronic interactions, under which the electrons in the system would reside
in a bound state and thus unable to propagate without overcoming the effective
binding energy.

Following Mott’s general picture, modern theoretical studies [12, 13] of Mott
insulators typically make use of half-filled Hubbard model [14, 15, 16, 17] or
its generalized form. For example, transition-metal oxides [6, 7, 8, 4] and
organic Mott insulators [18, 19] have relatively strong local electronic
repulsion $U$ compared to the kinetic energy. In half-filled systems, such a
large $U$ pushes states with doubly occupied orbitals toward high energy and
thus allow them to be “integrated out” from the low-energy sector, leaving
only effectively singly occupied states without charge freedom (thus an
insulator.) From the perspective of Mott’s picture, the repulsion $U$ in such
half-filled systems provides exactly the necessary electronic interaction to
bind the electron in an effective local orbital.

Still, an important physics question remains open, namely the role of the
long-range Coulomb interaction on the nature of the quantum Mott transition.
Mott argued [11, 1] that between the metallic and insulating phases the
significant change in the screening of long-range interaction could lead to a
first-order metal-insulator transition at zero temperature. In contrast,
current Hubbard model-based studies (which only incorporate local interaction)
found the quantum phase transition to be continuous [20, 12] instead. It is
therefore an important and timely task to investigate the role of long-range
interaction on the Mott transition, in particular to examine the validity of
Mott’s expectation.

Figure 1:  Mott’s picture of the key scales of metal-insulator transition: (a)
radius $a_{B}$ of the particle-hole bound state, and (b) range $\lambda$ of
screened long-range Coulomb attraction $V(r)$ as a function of relative
distance $r$. (c) When the bound state density $\rho_{B}$ is confined within
the Coulomb attraction, $a_{B}<\lambda$, the system is insulating. (d)
Otherwise the system becomes metallic. The finite $\lambda$ at the transition
dictates a first-order phase transition in this picture.

Here, to address this essential question, we attempt to incorporate the effect
of long-range interaction and investigate the Mott transition of a model
system with hydrogen atoms in a cubic lattice. The effect of interaction
screening is incorporated approximately via density self-consistency [21, 22]
of the density functional theory plus dynamic mean-filed theory (DFT+DMFT)
[23, 24, 25, 26, 27, 28, 29, 12, 13, 30]. We found that this model system is
in the charge-transfer regime, in which the relevant low-energy charge sector
is dominated by the 2$s$ and singly occupied 1$s$ orbitals. Correspondingly,
the Mott transition is intimately related to the closing of the charge-
transfer gap that indicates the destruction of possible bound states in the
metallic phase, in good consistency with Mott’s general argument. Importantly,
the insulating gap is found to reduce smoothly to zero and the charge carrier
density (reflected by $1s$ orbital occupations) grows continuously. Both
results suggest a continuous quantum phase transition, instead of the first-
order transition expected by Mott’s consideration of long-range screening. We
analytically shows that the key physics of this most common type of Mott
transition in materials involves only local correlation and local energy
scales, long-range screening therefore does not affect the nature of the
transition.

Figure 2:  (a) One-body density of states (DOS) in unit of (eV$\cdot$ unit
cell)-1 for $a=3.0\text{\AA}$ via standard DFT treatment, indicating a
metallic system with a charge-transfer gap $E_{CT}$ above the chemical
potential set as the reference energy at zero. (b) The same via DFT+DMFT at
$T=0.01~{}$eV, giving an insulating system with a correlation-renormalized
charge-transfer gap $\tilde{E}_{CT}$ across the chemical potential. (c) Smooth
reduction of charge-transfer gap $E_{CT}$ and $\tilde{E}_{CT}$ as the lattice
constant $a$ decreases toward the quantum critical point at $a_{c}$, below
which $\tilde{E}_{CT}<0$ and the system turns metallic with a gradual increase
of itinerant hole carrier density reflected by the smooth reduction of $1s$
occupation (d). Both (c) and (d) indicate a continuous quantum phase
transition different from Mott’s proposal.

Let’s first review the original consideration of the Mott transition. Mott
proposed that a lattice of one-electron atom, such as hydrogen, must be an
insulator at large lattice spacing. This is because starting from a ground
state with a uniform charge of one electron per atom, the low-energy current
response must be gapped associated with the formation of particle-hole bound
states with an effective Bohr radius $a_{B}$ [cf. Fig. 1(a)] (assuming a
poorly screened $e^{2}/r$ attractive Coulomb interaction over distance $r$.)
On the other hand, at small lattice spacing the system is expected to become a
good metal with a half-filled lowest energy band. This is because the current
response now turns gapless in association with the delocalization of particle-
hole pairs, as a result of the long-range screening that reduces the range of
the attraction $-(e^{2}/r)\text{exp}(-2\pi r/\lambda)$ to a finite length
scale $\lambda$ [cf. Fig. 1(b)]. Therefore, there must exist a transition
point for a metal-insulator transition at some intermediate lattice spacing.

As a simple estimation, Mott proposed that the metal-insulator transition
would take place at the point when $a_{B}$ and $\lambda$ are equal. This seems
quite reasonable since at $a_{B}\ll\lambda$ [cf. Fig. 1(c)] the bound state is
well preserved such that the system should be a good insulator, while at
$a_{B}\gg\lambda$ [cf. Fig. 1(d)] the interaction is screened too strongly to
maintain a bound state such that the system would be a good metal. Based on
this rough criteria, Mott estimated the critical density $n_{c}^{1/3}a_{B}\sim
0.25$ [31] via Thomas-Fermi approximation. Though this estimation is very
crude, the criteria actually works very well for many materials [2].

The necessity of a finite $\lambda\sim a_{B}$ in Mott’s criterion of metal-
insulator transition suggests strongly that the transition be first-order,
since a finite $\lambda$ corresponds to a finite itinerant carrier density. As
another way to visualize this, near the critical density, a slight increase of
itinerant electrons would lead to an enhanced screening (or reduced
$\lambda$), which in turn further increases the carrier density. One would
therefore expect that such a non-linear feedback can lead to a discontinuous
transition at zero temperature. Unfortunately due to the complexity of quantum
many-body problems, theoretical studies to date have not been able to
incorporate the long-range interactions, nor has the first-order quantum phase
transition been obtained [12, 13].

To include this important screening effect of long-range Coulomb interaction
beyond the Thomas-Fermi approximation, we use the density self-consistency
scheme [21, 22] of density functional theory (DFT) plus dynamical mean-field
theory (DMFT). In the process of self-consistent iteration, the change of
electron density and corresponding screening effect will feedback to the next
iteration and thus allow the possibility of the first-order transition. The
strong intra-atomic repulsion $U$ ($\sim 1$Ry [30]) is included via DMFT [25,
26, 27, 28, 29, 12, 13] which is so far the most widely applied approximation
to achieve the Mott insulating phase. This scheme is implemented by
DFT+embedded DMFT functional [21, 22].

We first start with a low-density insulating case at lattice spacing
$a=3.0~{}\text{\AA}$. Fig. 2(a) shows the resulting density of states (DOS) in
the standard DFT calculation via local density approximation [23, 24]. As
expected, the system is incorrectly identified as a metal with chemical
potential in the middle of the first band, corresponding to the half-filling
$1s$ orbital. Upon inclusion of strong local exchange-correlation via DMFT at
$T=0.01~{}$eV, Fig. 2(b) shows the correct insulating state with a DOS hosting
a gap around the chemical potential.

It is obvious that this insulating gap is not the standard Mott gap in a
Hubbard model, since its size [32] is only around $2~{}$eV in this case, much
smaller than the local repulsion $U\sim 1~{}$Ry. Indeed, in both DFT [panel
(a)] and DFT+DMFT [panel (b)] results, the orbital content right above the gap
is predominately of 2$s$ and 2$p$ instead of 1$s$. Therefore, within the
parameter range near the Mott transition, this model system is well within the
charge-transfer (CT) insulator regime, in which the lowest-energy charge
fluctuation is between 1$s$ and 2$s$/2$p$ orbitals of neighboring atoms.

Notice that the charge-transfer gap $E_{CT}$ is strongly reduced from $4~{}$eV
in the DFT result [panel(a)] to only $\sim 2~{}$eV in the DFT+DMFT result
[panel(b)]. This can be understood from the additional energy lowering from
the emergent antiferromagnetic spin exchange between the itinerant carrier in
2$s$/2$p$ and the local immobile 1$s$ spin, similar to the one that drives the
formation of Kondo singlet [33, 34].

Now approaching the metal-insulating transition via reducing the lattice
constant, as shown by the black line in Fig. 2(c), the CT gap of the
insulating solution reduces and eventually vanishes around $a=a_{c}\sim
2.8~{}\text{\AA}$. Naturally, after that the insulator solution can no longer
be found and the system is metallic. Incidentally, for $a>a_{c}$ metallic
solution is unstable as well. This indicates that the phase transition is
intimately related to the gap closing. (Note that in an extremely narrow range
of $\Delta a\sim 0.003~{}\text{\AA}$ around $a_{c}$, both metallic and
insulating solutions can be found [30]. Considerations on this nearly
negligible and somewhat artificial feature will be discussed later.)

Most importantly, our results display a very important deviation from Mott’s
expectation on the phase transition. Figure 2(c) shows that even with the
long-range Coulomb interaction and its screening incorporated via the DFT
treatment, as the lattice spacing decreases, $\tilde{E}_{CT}$ only smoothly
reduces to zero. Furthermore, as the lattice spacing decreases further, Fig.
2(d) shows that the occupation of the $1s$ orbital decreases smoothly as well,
indicating a gradual development of hole carrier density. No abrupt
enhancement of gap closing or the hole carrier density is found at $a_{c}$
that corroborates Mott’s expectation of a first-order phase transition.
Clearly, our results suggest a continuous metal-insulator transition. In other
words, within the approximation of the state-of-art DFT+DMFT approach
(omission of non-local correlation), the proposed drastic effect of screening
of long-range interaction does not appear to be essential to the Mott
transition in this model system.

One can analytically understand the weak relevance of long-range screening,
owing to the simplicity of the charge transfer nature of this model system.
First notice that in terms of the energy scale, the long-range screening
involves particle-hole attraction at long distances is at best on the order of
100 meV. In contrast, a realistic multi-orbital interacting Hamiltonian
obtained from Wannier function analysis [35, 36, 37, 38, 39] of the DFT result
indicates several stronger eV-scale physics. Therefore, before incorporating
effects of the long-range screening, one should consider first the following
representative dominant physics of our model system,

$\begin{split}H&=\tau\sum_{ii^{\prime}\mu}f^{\dagger}_{i\mu}f_{i^{\prime}\mu}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\\\
&+\sum_{jj^{\prime}\mu}t_{jj^{\prime}}c^{\dagger}_{j\mu}c_{j^{\prime}\mu}+\Delta\sum_{j\mu}c^{\dagger}_{j\mu}c_{j\mu}\\\
&+V\sum_{<ij>\mu}f^{\dagger}_{i\mu}c_{j\mu}+\text{H.c.},\end{split}$ (1)

where $f_{i\mu}^{\dagger}$ and $c^{\dagger}_{i\mu}$ are creation operators for
the 1$s$ and 2$s$ orbitals of spin $\mu$ located at site $i$, and
$n_{i\mu}=f^{\dagger}_{i\mu}f_{i\mu}$ is the density operator of 1$s$ orbital.
The first line describes kinetic hopping $\tau_{ii^{\prime}}$ of the heavily
occupied 1$s$ orbital and its strong intra-orbital repulsion $U$. The second
line states the kinetic hopping $t_{ii^{\prime}}$ of the slightly occupied
2$s$ orbital due to its higher site energy $\Delta$. The third line gives the
hybridization $V$ of 1$s$ orbital with the nearest neighboring 2$s$ orbital.
Naturally, the total electron count corresponds to one electron per hydrogen
atom,
$\sum_{i\mu}f^{\dagger}_{i\mu}f_{i^{\prime}\mu}+\sum_{j\mu}c^{\dagger}_{j\mu}c_{j\mu}=N$,
for a system with $N$ total sites.

As pointed out above, this model system turns out to be in the charge transfer
regime, $\tau<V\ll\Delta<U$. It is therefore natural to integrate out the
charge fluctuation that generates states with doubly occupied $1s$ orbital,
via a Schrieffer-Wolff transformation [30]. (This great simplification is not
possible in the standard Hubbard model, since its metallic states are all with
double occupation.) In the dilute carrier limit relevant to the quantum metal-
insulator transition of interest here, the low-energy effective Hamiltonian
thus can be written as:

$\begin{split}\tilde{H}&=J\sum_{<ii^{\prime}>}(\mathbf{S}_{i}\cdot\mathbf{S}_{i^{\prime}}-\frac{1}{4}\tilde{n}_{i}\tilde{n}_{i^{\prime}})+\sum_{ii^{\prime}\mu}\tilde{\tau}_{ii^{\prime}}\tilde{f}^{\dagger}_{i\mu}\tilde{f}_{i^{\prime}\mu}\\\
&+\sum_{jj^{\prime}\mu}\tilde{t}_{jj^{\prime}}\tilde{c}^{\dagger}_{j\mu}\tilde{c}_{j^{\prime}\mu}+\tilde{\Delta}\sum_{j\mu}\tilde{c}^{\dagger}_{j\mu}\tilde{c}_{j\mu}\\\
&+K\sum_{<ij>}(\mathbf{S}_{i}\cdot\sum_{\mu\nu}\tilde{c}^{\dagger}_{j\mu}\bm{\sigma}_{\mu\nu}\tilde{c}_{j\nu}-\frac{1}{2}\tilde{n}_{i}\sum_{\mu}\tilde{c}_{j\mu}^{\dagger}\tilde{c}_{j\mu})\\\
&+K\sum_{i,j\neq
j^{\prime},\mu}\mathbf{S}_{i}\cdot\sum_{\mu\nu}\tilde{c}_{j\mu}\bm{\sigma}_{\mu\nu}\tilde{c}_{j^{\prime}\nu}+L\sum_{i,j\neq
j^{\prime}\mu}\tilde{n}_{i}c_{j\mu}^{\dagger}c_{j^{\prime}\mu},\end{split}$
(2)

where $\tilde{n}_{i}=\sum_{\mu}\tilde{f}^{\dagger}_{i\mu}\tilde{f}_{i\mu}$ and
$\mathbf{S}_{i}=\frac{1}{2}\sum_{\mu\nu}\tilde{f}^{\dagger}_{i\mu}\bm{\sigma}_{\mu\nu}\tilde{f}_{i\nu}$
are the density and spin operator of the low-energy renormalized particle
$\tilde{f}^{\dagger}_{i\mu}$. $\bm{\sigma}_{\mu\nu}$ are the Pauli matrices.

For the low-energy metallic states containing dilute carriers, the heavily
renormalized kinetic strength $\tilde{\tau}$ of the 1$s$ orbitals,
$\tilde{f}^{\dagger}_{i}$, is only active for the dilute hole carriers
residing in them, since they no longer allow double occupation, i.e.
$\tilde{f}^{\dagger}_{i\uparrow}\tilde{f}^{\dagger}_{i\downarrow}=0$.
Correspondingly, most 1$s$ electrons behave like immobile local magnetic
moments coupling to each other via $J=4\tau^{2}/U$. The corresponding strong
magnetic correlation in turn constrains the hole motion as in the $t$-$J$
model. In contrast, the kinetic strength $\tilde{t}$ of the dilute electron
carriers residing in the renormalized 2$s$ orbitals,
$\tilde{c}^{\dagger}_{j}$, is weakly renormalized, but they develop a
correlated hopping $K=V^{2}(\frac{1}{U-\Delta}+\frac{1}{\Delta})$
[$L=\frac{1}{2}V^{2}(\frac{1}{\Delta}-\frac{1}{U-\Delta})$] sensitive to the
spin direction of the 1$s$ local moments, which would promote Kondo-like
physics such as the Kondo resonant peak at the chemical potential. On the
other hand, for insulating states containing zero electron carriers
$\tilde{c}^{\dagger}_{j}\tilde{c}_{j}=0$ and
$\tilde{f}^{\dagger}_{i}\tilde{f}_{i}=1$, only the first term of Eq. 2 is
active and the effective Hamiltonian reduces to the standard Heisenberg model.

Figure 3: Comparison of the insulating mechanism in Mott’s picture and our
model system. (a) In Mott’s picture, particle and hole carriers are assumed to
be bound by Coulomb attraction when the length scale of long-range screening
$\lambda$ is larger than the size of the bound state. (b) In our model system,
originating from the strong intra-atomic repulsion $U$, the effective
attraction for particles in the 1$s$ orbitals is only within the length scale
of lattice spacing $a$. The metal-insulator transition involves instead
particles in the 2$s$ orbitals, corresponding to the closing of the charge-
transfer gap $\tilde{\Delta}-W/2$, with $W$ being the sum of bandwidth of
dressed particles and holes. These dominant factors are all local in space,
the long-range screening does not play an essential role.

Now it is clear why the long-range screening plays little role in this model
system (and in fact most of the real Mott insulating materials). In charge-
transfer systems, $\tau<V\ll\Delta<U$, the insulating states are stabilized by
the strong intra-orbital repulsion $U$ and orbital energy difference $\Delta$.
Specifically, one can create a particle-hole pair residing within the 1$s$
orbitals in different atoms, but the electron carriers would necessarily lead
to a double occupation. This gives an effective short-range particle-hole
binding potential $V_{B}\sim-U\theta(\lambda-r)$, where $\lambda$ is just the
atomic radius shown as the orange line in Fig. 3(b). Since the binding is
dominated by the short-range physics, long-range screening naturally bears
little relevance in this case.

Alternatively, one can create a lower-energy particle-hole (2$s$-1$s$) pair,
which requires at least the charge transfer energy
$\sim\tilde{\Delta}-6(\tilde{t}+\tilde{\tau})=\tilde{\Delta}-W/2$ ($W$ being
the sum of bandwidth of particle and hole bands, subject to all possible
further renormalization at lower-energy scale.) [cf. the black dash line in
Fig. 3(b)]. Upon reducing the inter-atomic distance $a$ (and increasing
$\tilde{t}$ and $\tilde{\tau}$), this energy would eventually turn negative
and the ground state would host a small number of particles and an equal
amount of holes to further lower the energy.

Note, however, that these particle-holes can still bind weakly as charge-
neutral excitons such that the system can still remain insulating as an
excitonic insulator [40, 41] until further enhancement of the bare kinetic
energy demands a larger number of particles and holes with Fermi energies
exceeding the weak particle-hole binding. At that point, the system will
finally become truly metallic.

In such an excitonic insulating phase, since the exciton binding energy is
typically on the order of 100 meV or lower, one might still consider the
relevance of Mott’s argument of long-range screening. Nonetheless, notice that
in our model system, these two types of carriers reside on different atomic
orbitals of very different sizes, and their kinetics are renormalized by
distinct mechanisms discussed above. It is therefore unlikely for the particle
$\tilde{c}^{\dagger}$ and hole $\tilde{f}^{\dagger}$ carrier, having very
different velocity and effective mass, to bind in the typical energy scale.
Consequently, such an excitonic insulating phase, in which the long-range
screening is able to play some role, must occur in a negligibly narrow region,
if even exists, such that the MIT should take place in close proximity to the
vanishing of the charge transfer energy as found in Fig. 2(c).

Finally, it is well-known [12, 13, 42] that the lack of inter-atomic
correlation in the DMFT treatment can lead to unphysically large degree of
degeneracy and consequently artificially large entropy at finite temperature
$T=0^{+}$. This often causes an entropy-driven first-order transition at
finite temperature [43, 44]. Our finite-temperature calculation and the
effective analytical model above are not immune to this problem [30].
Nonetheless, it is interesting to note that this artifact is much weaker in
the charge-transfer systems than that in the Hubbard model. In fact, in Fig.
2(c), the region hosting first-order transition is narrower than the thickness
of the black line. We take it that this distinction reflects the more robust
continuous nature of the quantum phase transition in such charge-transfer
systems. Also, although unlikely from an energy scale point of view, future
investigation on potential subtle roles of long-range interaction and its
screening via more thorough treatments such as DFT+$GW$+EDMFT [45, 46, 47]
might still be of interest.

Experimentally, our model system of periodic hydrogen lattice could be tested
by a new generation of capabilities that engineers hydrogen-like lattice of
any desired geometry through STM-implanting donors on surfaces of
semiconductors [48]. These experimental advances allow tuning the system
across Mott transition without randomness and are thus ideal to experimentally
clarify the role of long-range interaction.

In summary, we investigate the important open question concerning the role of
long-range screening on the quantum Mott transition using a model system of
hydrogen lattice. We employ the charge self-consistent scheme of the dynamical
mean-field theory that incorporates approximately the long-range interaction
and its screening within the density functional treatment. The system is found
to be well within the charge-transfer regime and the charge-transfer gap
closes smoothly across the quantum phase transition. Together with a gradually
increasing hole carrier density in the metallic phase, this indicates a
continuous quantum phase transition distinct from Mott’s proposal. We explain
this deviation from Mott’s long-range screening-based first-order picture by
demonstrating analytically that 1) the obtained insulating phase in this
typical case is driven by short-range physics, and 2) the metal-insulator
transition involves closing of the charge transfer gap controlled by only
local energy scales. Our result suggests that the metal-insulator transitions
in a large class of real-life materials are not governed by Mott’s (sufficient
but not necessary) long-range screening-based scenario.

We thank Yuting Tan and Xinyi Li for their helpful discussions. This work is
supported by National Natural Science Foundation of China (NSFC) under Grants
No. 12274287 & 12042507, and Innovation Program for Quantum Science and
Technology No. 2021ZD0301900. S.S. acknowledges support from the Science and
Engineering Research Board, India (SRG/2022/000495), (MTR/2022/000638), and
IIT(ISM) Dhanbad [FRS(175)/2022- 2023/PHYSICS]. Work in Florida (P. K. H. T.
and V. D.) was supported by the NSF Grant No. 1822258, and the National High
Magnetic Field Laboratory through the NSF Cooperative Agreement No. 1644779
and the State of Florida.

## References

  * MOTT [1968] N. F. MOTT, Metal-insulator transition, Rev. Mod. Phys. 40, 677 (1968).
  * Mott [1990] N. Mott,  _Metal-Insulator Transitions_, 1st ed. (CRC Press, London, 1990).
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Dagotto [1994] E. Dagotto, Correlated electrons in high-temperature superconductors, Rev. Mod. Phys. 66, 763 (1994).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * de Boer and Verwey [1937] J. H. de Boer and E. J. W. Verwey, Semi-conductors with partially and with completely filled 3d-lattice bands, Proceedings of the Physical Society 49, 59 (1937).
  * Brito  _et al._ [2016] W. H. Brito, M. C. O. Aguiar, K. Haule, and G. Kotliar, Metal-insulator transition in ${\mathrm{vo}}_{2}$: A $\mathrm{DFT}+\mathrm{DMFT}$ perspective, Phys. Rev. Lett. 117, 056402 (2016).
  * Zhang  _et al._ [2017] H. Zhang, K. Haule, and D. Vanderbilt, Metal-insulator transition and topological properties of pyrochlore iridates, Phys. Rev. Lett. 118, 026404 (2017).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, Rev. Mod. Phys. 75, 473 (2003).
  * Anderson [1950] P. W. Anderson, Antiferromagnetism. theory of superexchange interaction, Phys. Rev. 79, 350 (1950).
  * Mott [1949] N. F. Mott, The basis of the electron theory of metals, with special reference to the transition metals, Proceedings of the Physical Society. Section A 62, 416 (1949).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
  * Hubbard and Flowers [1963] J. Hubbard and B. H. Flowers, Electron correlations in narrow energy bands, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Hubbard and Flowers [1964a] J. Hubbard and B. H. Flowers, Electron correlations in narrow energy bands. II. The degenerate band case, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 277, 237 (1964a).
  * Hubbard and Flowers [1964b] J. Hubbard and B. H. Flowers, Electron correlations in narrow energy bands III. An improved solution, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 281, 401 (1964b).
  * Gutzwiller [1963] M. C. Gutzwiller, Effect of correlation on the ferromagnetism of transition metals, Phys. Rev. Lett. 10, 159 (1963).
  * Powell and McKenzie [2011] B. J. Powell and R. H. McKenzie, Quantum frustration in organic mott insulators: from spin liquids to unconventional superconductors, Reports on Progress in Physics 74, 056501 (2011).
  * Pustogow  _et al._ [2018] A. Pustogow, M. Bories, A. Löhle, R. Rösslhuber, E. Zhukova, B. Gorshunov, S. Tomić, J. A. Schlueter, R. Hübner, T. Hiramatsu, Y. Yoshida, G. Saito, R. Kato, T. H. Lee, V. Dobrosavljević, S. Fratini, and M. Dressel, Quantum spin liquids unveil the genuine mott state, Nature Materials 17, 773 (2018).
  * Rozenberg  _et al._ [1994a] M. J. Rozenberg, G. Moeller, and G. Kotliar, The metal-insulator transition in the hubbard model at zero temperature ii, Modern Physics Letters B 08, 535 (1994a).
  * Haule  _et al._ [2010] K. Haule, C.-H. Yee, and K. Kim, Dynamical mean-field theory within the full-potential methods: Electronic structure of ${\text{CeIrIn}}_{5}$, ${\text{CeCoIn}}_{5}$, and ${\text{CeRhIn}}_{5}$, Phys. Rev. B 81, 195107 (2010).
  * Haule and Birol [2015] K. Haule and T. Birol, Free energy from stationary implementation of the $\mathrm{DFT}+\mathrm{DMFT}$ functional, Phys. Rev. Lett. 115, 256402 (2015).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Correlated lattice fermions in $d=\infty$ dimensions, Phys. Rev. Lett. 62, 324 (1989).
  * Müller-Hartmann [1989] E. Müller-Hartmann, Fermions on a lattice in high dimensions, International Journal of Modern Physics B 03, 2169 (1989).
  * Brandt and Mielsch [1989] U. Brandt and C. Mielsch, Thermodynamics and correlation functions of the falicov-kimball model in large dimensions, Zeitschrift für Physik B Condensed Matter 75, 365 (1989).
  * Janiš [1991] V. Janiš, A new construction of thermodynamic mean-field theories of itinerant fermions: application to the falicov-kimball model, Zeitschrift für Physik B Condensed Matter 83, 227 (1991).
  * Georges and Kotliar [1992] A. Georges and G. Kotliar, Hubbard model in infinite dimensions, Phys. Rev. B 45, 6479 (1992).
  * [30] See supplemental materials at xxx for all calculation details and discussions on the phase diagram.
  * Mott [1956] N. F. Mott, On the Transition to Metallic Conduction in Semiconductors, Canadian Journal of Physics 34, 1356 (1956).
  * [32] See supplemental materials at xxx for determination of charge-transfer gap $E_{CT}$ and renormalized charge-transfer gap $\tilde{E}_{CT}$.
  * Kondo [1964] J. Kondo, Resistance Minimum in Dilute Magnetic Alloys, Progress of Theoretical Physics 32, 37 (1964).
  * Ng and Lee [1988] T. K. Ng and P. A. Lee, On-site coulomb repulsion and resonant tunneling, Phys. Rev. Lett. 61, 1768 (1988).
  * Lang  _et al._ [2021] Z.-J. Lang, R. Jiang, and W. Ku, Strongly correlated doped hole carriers in the superconducting nickelates: Their location, local many-body state, and low-energy effective hamiltonian, Phys. Rev. B 103, L180502 (2021).
  * Jiang  _et al._ [2022] R. Jiang, Z.-J. Lang, T. Berlijn, and W. Ku, Variation of carrier density in semimetals via short-range correlation: a case study with nickelate NdNiO2, arXiv e-prints , arXiv:2207.08802 (2022), arXiv:2207.08802 .
  * Lee  _et al._ [2013] C.-C. Lee, X. M. Chen, Y. Gan, C.-L. Yeh, H. C. Hsueh, P. Abbamonte, and W. Ku, First-principles method of propagation of tightly bound excitons: Verifying the exciton band structure of lif with inelastic x-ray scattering, Phys. Rev. Lett. 111, 157401 (2013).
  * Yin  _et al._ [2006] W.-G. Yin, D. Volja, and W. Ku, Orbital ordering in ${\mathrm{lamno}}_{3}$: Electron-electron versus electron-lattice interactions, Phys. Rev. Lett. 96, 116405 (2006).
  * Ku  _et al._ [2002] W. Ku, H. Rosner, W. E. Pickett, and R. T. Scalettar, Insulating Ferromagnetism in ${\mathrm{L}\mathrm{a}}_{4}{\mathrm{B}\mathrm{a}}_{2}{\mathrm{C}\mathrm{u}}_{2}{\mathrm{O}}_{10}$: An Ab Initio Wannier Function Analysis, Phys. Rev. Lett. 89, 167204 (2002).
  * Jérome  _et al._ [1967] D. Jérome, T. M. Rice, and W. Kohn, Excitonic insulator, Phys. Rev. 158, 462 (1967).
  * Kohn [1967] W. Kohn, Excitonic phases, Phys. Rev. Lett. 19, 439 (1967).
  * Vollhardt  _et al._ [2012] D. Vollhardt, K. Byczuk, and M. Kollar, Dynamical mean-field theory, in  _Strongly Correlated Systems: Theoretical Methods_, edited by A. Avella and F. Mancini (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 203–236.
  * Georges and Krauth [1993] A. Georges and W. Krauth, Physical properties of the half-filled hubbard model in infinite dimensions, Phys. Rev. B 48, 7167 (1993).
  * Rozenberg  _et al._ [1994b] M. J. Rozenberg, G. Kotliar, and X. Y. Zhang, Mott-hubbard transition in infinite dimensions. ii, Phys. Rev. B 49, 10181 (1994b).
  * Biermann  _et al._ [2003] S. Biermann, F. Aryasetiawan, and A. Georges, First-Principles Approach to the Electronic Structure of Strongly Correlated Systems: Combining the $GW$ Approximation and Dynamical Mean-Field Theory, Phys. Rev. Lett. 90, 086402 (2003).
  * Ayral  _et al._ [2013] T. Ayral, S. Biermann, and P. Werner, Screening and nonlocal correlations in the extended Hubbard model from self-consistent combined $GW$ and dynamical mean field theory, Phys. Rev. B 87, 125149 (2013).
  * Petocchi  _et al._ [2020] F. Petocchi, F. Nilsson, F. Aryasetiawan, and P. Werner, Screening from ${e}_{g}$ states and antiferromagnetic correlations in ${d}^{(1,2,3)}$ perovskites: A $GW+\text{EDMFT}$ investigation, Phys. Rev. Research 2, 013191 (2020).
  * Schofield  _et al._ [2013] S. R. Schofield, P. Studer, C. F. Hirjibehedin, N. J. Curson, G. Aeppli, and D. R. Bowler, Quantum engineering at the silicon surface using dangling bonds, Nature Communications 4, 1649 EP (2013).