11institutetext: Leiden Observatory, Leiden University, PO Box 9513, 2300 RA
Leiden, The Netherlands. 11email<EMAIL_ADDRESS>22institutetext:
Instituto de Radioastronomía y Astrofísica (IRyA), Universidad Nacional
Autónoma de México (UNAM). Campus Morelia (Michoacán, México).
33institutetext: ESO Garching, Karl-Schwarzschild-Str. 2, 85748, Garching bei
Munchen, Germany. 44institutetext: Department of Physics and Astronomy,
University of Exeter, Exeter, EX4 4QL, UK. 55institutetext: Mullard Space
Science Laboratory, University College London, 3 Holmbury St Mary, Dorking,
Surrey RH5 6NT, UK. 66institutetext: Jet Propulsion Laboratory, California
Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
77institutetext: Dipartimento di Fisica, Università degli Studi di Milano, Via
Giovanni Celoria 16, 20133 Milano, Italy 88institutetext: Joint ALMA
Observatory, Alonso de Córdova 3107, Vitacura, Santiago, 763 0355, Chile.
99institutetext: National Astronomical Observatory of Japan, Los Abedules 3085
Oficina 701, Vitacura, Santiago, 763 0000, Chile. 1010institutetext: The
Graduate University for Advanced Studies, SOKENDAI, Osawa 2-21-1, Mitaka,
Tokyo 181-8588, Japan.

# Into the thick of it:
ALMA 0.45 mm observations of HL Tau at 2 au resolution

Osmar M. Guerra-Alvarado  Carlos Carrasco-González 11 Enrique Macías 22 Nienke
van der Marel 33 Adrien Houge 11 Luke T. Maud 3344 Paola Pinilla 33 Marion
Villenave 55 Yoshiharu Asaki Elizabeth Humphreys 66778899101033

###### Abstract

Aims. To comprehend the efficiency of dust evolution within protoplanetary
disks, it is crucial to conduct studies of these disks using high-resolution
observations at multiple wavelengths with the Atacama Large
Millimeter/submillimeter Array (ALMA).

Methods.  In this work, we present high-frequency ALMA observations of the HL
Tau disk using its Band 9 centered at a wavelength of 0.45 mm. These
observations achieve the highest angular resolution in a protoplanetary disk
to date, 12 milliarcseconds (mas), allowing the study of the dust emission at
scales of 2 au. We use these data to extend the previously published multi-
wavelength analysis of the HL Tau disk, constraining the dust temperature,
dust surface density, and maximum grain size throughout the disk. We perform
this modeling for compact solid dust particles, as well as for porous
particles.

Results. Our new 0.45 mm data traces mainly optically thick emission,
providing a tight constraint to the dust temperature profile. We derive
maximum particle sizes of $\sim 1$ cm from the inner disk to $\sim 60$ au.
Beyond this radius, we find particles between 300 $\mu$m and 1 mm. The total
dust mass of the disk is 2.1 MJ with compact grains, and it increases to 6.3
MJ assuming porous particles. Moreover, an intriguing asymmetry is observed at
32 au in the northeast inner part of the HL Tau disk at 0.45 mm. We propose
that this asymmetry is the outcome of a combination of factors including the
optically thick nature of the emission, the orientation of the disk, and a
relatively large dust scale height of the grains preferentially traced at 0.45
mm. To validate this, we conducted a series of radiative transfer models using
the RADMC-3D software. Our models, incorporating varying dust masses and scale
heights, successfully replicate the observed asymmetry in the HL Tau disk. If
this scenario is correct, our measured dust mass within 32 au would suggest a
dust scale height H/R$>0.08$ for the inner disk. Finally, the unprecedented
resolution allowed us to probe for the first time the dust emission down to a
few au scales. We observed an increase in brightness temperature inside the
estimated water snowline and speculate whether it could indicate the presence
of a traffic jam effect in the inner disk.

Conclusions. Our results showcase how 0.45 mm observations of protoplanetary
disks can be used to robustly constrain their dust temperature radial profile.
Additionally, the higher optical depths at this wavelength can be used to
constrain the dust vertical scale height. Finally, these higher frequencies
allow us to reach higher spatial resolutions that have the potential to
resolve the region within the water snowline in disks.

###### Key Words.:

Planetary systems: Protoplanetary disk - Radio continuum: planetary systems

##  1 Introduction

Over the past ten years, the Atacama Large Millimeter/submillimeter Array
(ALMA) has made significant contributions to our understanding of
protoplanetary disks by revealing substructures within them. These studies
(e.g. ALMA Partnership et al. 2015, Isella et al. 2016, and Long et al. 2018)
have shown the frequent appearance of substructures in protoplanetary disks
having different morphologies (van der Marel et al. 2013, Pérez et al. 2016,
and Huang et al. 2018). These substructures play a significant role in
protoplanetary disks allowing dust grains to accumulate and grow, facilitating
the formation of planets by providing regions of higher dust density and
trapping, and slowing down the radial drift of dust grains, which would
otherwise hinder the planet formation process (Pinilla et al. 2012a).

Most ALMA studies on protoplanetary disks are performed at wavelengths around
1 mm, where dust emission is relatively bright and high angular resolutions
around 30 mas can be achieved. In some disks, similar studies have also been
performed at longer wavelengths where the emission is known to be optically
thinner, up to 3 mm with ALMA (e.g., Tazzari et al. 2021, Cazzoletti et al.
2018,Long et al. 2020) but also up to 1 cm with the Very Large Array (VLA)
(e.g., Carrasco-González et al. 2016, Hashimoto et al. 2023). These
observations at long wavelengths allow us to study the properties of the dust
grains albeit at the expense of angular resolution. Observations of
protoplanetary disks with ALMA at shorter wavelengths, although providing the
highest angular resolutions possible with ALMA ($\sim$10 mas), are usually
avoided because of the higher atmospheric opacity and turbulence, leading to
necessarily longer integration and calibration times, and the fact that dust
emission is expected to be optically thick. However, there is also
complementary information that can be extracted from studies at such
wavelengths. For instance, optically thick emission impose better constraints
on the temperature and the albedo of the dust particles, improving the results
from multi-wavelength modeling of disks. Moreover, observations at short
wavelengths are the only way to explore the dust content in the gaps whose
emission is usually very optically thin to be detected at long wavelengths.

Recent multi-wavelength analyses of Class II protoplanetary disks (Carrasco-
González et al. 2019, Macías et al. 2021, Sierra et al. 2021, and Guidi et al.
2022) have measured dust properties within these substructures. The dust sizes
determined from these observations are around 1 mm. However, there is a
discrepancy with measurements obtained from polarization observations (Kataoka
et al. 2017) which suggest dust sizes on the order of a few hundred
micrometers. To reconcile these conflicting results, some studies by Kataoka
et al. 2017 and Ohashi et al. 2020 have proposed the existence of two dust
populations. The larger particles (mm-sized) would experience more efficient
vertical settling, becoming decoupled from the gas and gathering towards the
midplane of the disk. On the other hand, smaller particles would remain mixed
with the gas in higher layers of the disk (Barrière-Fouchet et al. 2005). More
recently, Zhang et al. 2023 proposed that large (sizes larger than 1 mm) and
porous (porosities higher than 70$\%$) particles could reproduce both, multi-
wavelength and polarization results. A very recent polarization
multiwavelength study strongly supported the presence of porous particles in
HL Tau (Lin et al., 2023a)

Studies by Pinte et al. 2016, Villenave et al. 2020, Doi & Kataoka 2021,
Villenave et al. 2022 and Pizzati et al. 2023 have focused on understanding
the vertical structure of protoplanetary disks. Although limited in number,
these studies suggest that Class II disks exhibit significantly low dust scale
heights, particularly at longer wavelengths observed with ALMA. This is a
consequence of the distribution of dust grains changing along the vertical
direction due to settling (Dullemond & Dominik 2004,Pinilla et al. 2021). Most
constraints presented in these papers remain as upper limits (except in
HD163296), with suggestions such as H(100au) ¡ 1 au. Additionally, Pizzati et
al. 2023 discovered values ¡ 4 au for certain DSHARP disks, particularly those
with less favorable orientations. These findings were predominantly derived
from the presence/depth of gaps and rings in the outer regions of the disks
(outer 100 au) where the scale height in the inner regions remains mostly
unknown.

In this work, we present high-resolution Band 9 observations of HL Tau, a
young stellar object in the Taurus-Auriga molecular cloud ($\sim 147$pc, Galli
et al. 2018) with an estimated age of less than 1 million years (Liu et al.
2017). HL Tau is considered a Class I-II disk (Furlan et al. 2008) and it was
the first protoplanetary disk in which the presence of dark and bright rings
was discovered (ALMA Partnership et al. 2015). Since it has been intensively
observed in the last years with ALMA and VLA, high-quality data covering a
wide range of wavelengths, between 0.9 mm and 1 cm, are currently available
(Carrasco-González et al. 2016, Carrasco-González et al. 2019). There have
been also several previous studies modeling its continuum emission to infer
the dust properties (e.g., Jin et al. 2016, Pinte et al. 2016, Liu et al.
2017). From these, it has been well established that the dark rings are
actually gaps in the dust distribution, i.e. regions of lower density than the
adjacent bright rings. The origin of the substructures is still debated, but
the most promising are planet-disk interactions (e.g., Jin et al. 2016,
Dipierro et al. 2016). Our new sensitive and high-resolution observations at
0.45 mm data are combined with the previous multi-wavelength observations to
analyze the dust properties in the disk. The addition of shorter wavelength
observations allowed us to improve previous modeling of the dust properties in
the disk, as well as obtain new insights about the dust temperature very close
to the protostar and the vertical structure of the disk.

Figure 1:  HL Tau Band 9 images. Left panel: Image using a robust parameter of
0 with a beam of 0.0226”$\times$ 0.021” Right-top panel: Zoom into the inner
part of the HL Tau disk of the robust 0 images showing the asymmetry with gray
contours at 11 times the rms and 33 times the rms. Right-bottom panel: Inner
HL Tau disk from the robust -0.5 image with a beam of 0.0116”$\times$ 0.010”.
Figure 2: Comparison of the HL Tau Band 6-7 (ALMA Partnership et al. 2015) and
the HL Tau Band 9 radial profiles of the brightness temperature. In yellow we
plotted the Band 6-7 ALMA image with the lowest resolution (33 mas), In black
the Band 9 image with 22 mas, and in green the Band 9 image with 13 mas. The
substructures are shallower in the radial brightness temperature for Band 9
compared to the Band 6-7 image, moreover, the emission from the Band 9 image
extends up to 140 au. In the inside panel, we show the first 20 au of the HL
Tau disk, we plotted the brightness temperatures of the three images and,
assuming that Tdisk=TB because of the optically thick regime, we displayed the
position (in temperature) of the water snow line.

##  2 Observations

We analyze archival observations obtained with the Atacama Large
Millimeter/submillimeter Array (ALMA) at a wavelength of 0.45 mm. One of the
archival datasets that we used was part of the test observations for the band-
to-band calibration during the High-Frequency Long Baseline Campaign (Asaki et
al., 2020), project code: 2011.0.00005.E). The observations were taken during
two execution blocks on November 3, 2017. The total observing time was 65 and
100 minutes for the first and second execution, respectively, which combined
gives a 45-minute on-source time for HL Tau. The final data has 8 spectral
windows, each of them covering a bandwidth of 2 GHz. We also used additional
data at a lower resolution (project code: 2017.1.01178.S, PI: E. Humphreys) to
complete the uv coverage at shorter baselines. These data have 8 spectral
windows covering frequencies from 646-662 GHz with a bandwidth of 1.875 GHz
each. The observations were taken from October 2, 2018 to October 18, 2018
with a total observing time on-source of 95 minutes.

We calibrated the data by using the calibration pipeline and scripts provided
by the ALMA staff. We used version 5.6.1 of CASA (Common Astronomy Software
Applications; McMullin et al. 2007 for self-calibration and cleaning in both
data sets. We identified and flagged the spectral lines that were found in
both data sets. We also averaged channels for the low and high-resolution
data; the final self-calibrated datasets contain spectral windows of 8
channels of 200 MHz each.

We performed phase-only and amplitude self-calibration in both data sets.
First, self-calibration was performed on the short baseline data. During the
self-calibration process, we created the images using the task clean in CASA.
The MTMFS deconvolver (Rau & Cornwell, 2011) was used with scales of 0, 1, 3,
and 5 times the beam size. We also used a Briggs weighting with a robust
parameter of 0.5. For the low-resolution data set the self-calibration process
resulted in a big improvement in the image quality. We applied 6 iterations of
phase-only self-calibration decreasing the time solution interval down to 18s.
Furthermore, we were also able to apply two rounds of amplitude and phase
self-calibration (calmode=ap) until no improvement in the signal-to-noise
ratio was found. By the end of the whole self-calibration process, we achieved
an increase in the S/N by a factor of $\sim$17.4. We combined our short
baseline data together with the long baseline data after aligning and
correcting for proper motions and astrometric errors between our data sets.
This was done by applying a Gaussian fit using CASA in the high-resolution
image and using the tasks fixvis and fixplanets to apply shifts to the low-
resolution data. After this, self-calibration was attempted on the combined
dataset. We started with phase-only self-calibration, however, due to the low
S/N at the long baselines, no improvement was found.

Lastly, self-calibration separately on the long-baseline data was also
attempted by using the model obtained from the self-calibrated low-resolution
data, but no improvement was found either. Since not enough solutions were
found with any method to improve the quality of the combined data set, the
self-calibration process was stopped.

With our final data set at hand, we recalculated the weights using the task
statwt in CASA. For the final images, we used the tclean task in CASA with the
MTMFS deconvolver (Rau & Cornwell, 2011) using nterms=1 and scales at 0, 1, 3,
and 5 times the beam size. We used a pixel size 10 times smaller than the beam
size for all of the images. We explored different weightings: natural,
uniform, and Briggs with different robust values. Finally, two images are used
for this study. The first one, made by setting the parameter robust to 0 and
using baselines shorter than 12 $k\lambda$ in length, gives us a good
compromise between resolution and sensitivity and was used for the study of
the large scale emission of the disk. This image has an rms noise of 270
$\mu$Jy beam-1 and a beam size of $\sim$23 mas. A second image was done by
setting the parameter robust to -0.5, which results in a beam size of $\sim$12
mas and an rms noise of 380 $\mu$Jy beam-1. This higher angular resolution
image allowed us a detailed study of the innermost part of the disk
($\lesssim$20 au). Both images are shown in Figure 1.

In order to study the radial properties of the disk, we also obtained radial
profiles of the emission and convolved the images to a round beam (robust 0,
24 mas and robust -0.5, 13 mas). For this, we used the known inclination and
position angles of the disk, 46.72∘ and 138.02∘, respectively (ALMA
Partnership et al., 2015). However, due to the inclination of the disk, there
is a loss in angular resolution in the NE-SW direction. Since the S/N of the
image is very high, we could avoid this loss in resolution by simply averaging
emission azimuthally within $\pm$0.2 rad of the major axis of the disk. Radial
profiles of the brightness temperature obtained from the intensity profiles of
both images are shown in Figure 2.

##  3 Results

In Figures 1 and 2, it can be seen that our images at 0.45 mm recover emission
from the whole disk with a radius of $\sim$1$\arcsec$ or $\sim$140 au, similar
to previous ALMA images at longer wavelengths (see e.g., Carrasco-González et
al. 2019, ALMA Partnership et al. 2015). Our 0.45 mm image also shows some
substructures but, as expected, due to the higher optical depth at this
wavelength, substructures are detected with a lower contrast between dark and
bright rings. This suggests that the gaps contain small particles traced by
the 0.45 mm image, these particles are more coupled with the gas and might
fill the gaps more than bigger particles observed at longer wavelengths
(Pinilla et al. 2012b,de Juan Ovelar et al. 2013). Furthermore, we observed a
slight deviation in the positions of the gaps and rings at 0.45 mm when
compared to the radial profiles at other wavelengths (see Figure 2), which
suggests that there might be slight variations in the depth and width of the
disk substructures at different wavelengths.

In the right-top of Figure 1, we show a zoomed-in view of the inner 50 au. An
asymmetry can be seen in the 0.45 mm data that is not seen at longer
wavelengths: the NE part of the first ring is more intense than the SW part of
the ring.

In the right-bottom panel of Figure 1, we show our highest angular resolution
image which allowed us to resolve, for the first time, the most internal part
of the disk ($\lesssim 20$ au). In this image, we can observe the first ring
and a very deep gap.

This first gap, centered at 13 au, while it shows some emission in the low
resolution 0.45 mm, it appears very dark in our highest angular resolution
image (see Figures 1 and 2).

These observations provide a clearer view of the gap’s structure, revealing
that it is likely deeper and more depleted of dust particles than what was
initially observed with lower resolution at 0.87 mm and 1.3 mm. Moreover, the
fact that we also detect a very deep gap in the 0.45 mm image strongly
suggests that this gap is actually devoid of small dust particles as well.

In the inner panel of Figure 2, it can be seen that the brightness temperature
increases within 2.5 au when observed at the highest angular resolution. This
could potentially be attributed to substructures within the inner disk that
are becoming resolved for the first time thanks to the higher angular
resolution of our 0.45 mm data.

In the following section, we first discuss the updated modeling of the radial
dust properties using all the currently available data of HL Tau, then the
origin of the new features found in the 0.45 mm images, the asymmetry, and the
substructure within the inner disk.

##  4 Discussion

###  4.1 Modeling of the dust properties

Figure 3: Dust parameters at each radius using compact dust particle
opacities, the position of the most common snowlines are plotted in colors.
The vertical dotted and dashed lines are the positions of the gaps and rings,
respectively, as seen from the surface density parameter. Top panel: Dust
temperature at each radius, the temperature is very well constrained due to
the optically thick nature of the Band 9 image. Middle panel: Dust surface
density at each radius. As expected the density increases in the bright rings
and decreases in the gaps. Bottom panel: Maximum particle size at each radius,
the particle sizes are around 1 mm - 1 cm.  Figure 4: Same as Figure 5 but
using a 90% porosity for the dust particles and for the dust opacities. The
dust temperature is lower, the surface density is higher around the gaps and
the maximum particle sizes are similar, however, at the 13 au gap there is a
huge increase in the particle size.

Previous to this work, there have been several studies modeling the dust
properties in the HL Tau using the available data at that moment (e.g., Pinte
et al. 2016, Liu et al. 2017, Carrasco-González et al. 2019). We perform a new
modeling procedure by including the new 0.45 mm image reported here which will
provide valuable constraints on the dust temperature and will enhance the
comprehensiveness of our multi-wavelength analysis. Thus, this new analysis
includes images at 0.45, 0.87, 1.3, 2.1 and 7.0 mm. We created a new 0.45 mm
image with a robust parameter of 1.0, using baselines ¡ 12 $k\lambda$ in
length as previously described for the robust 0.0 image. Then, we convolved it
with a 50 mas Gaussian beam and extracted radial profiles in wedges along the
major axis within radial beams of approximately 0.01 arcseconds. We also
extracted radial profiles at the other wavelengths applying the same procedure
to the images presented in Carrasco-González et al. (2019). Subsequently, each
point in the radial profiles was modeled independently. We note that each
point in the model is not independent from its neighboring points due to the
effects of the beam. Consequently, the contrast of the model parameters
between rings and gaps will be slightly smoothed in our results. We follow a
similar approach to (Carrasco-González et al., 2019), i.e., at each radius, we
fit the spectral energy distribution (SED) in the millimeter wavelength range
using a 1D slab model given by:

$\displaystyle
I_{\nu}=B_{\nu}(T)[(1-\exp(\tau_{\nu}/\mu))+\omega_{\nu}F(\tau_{\nu},\omega_{\nu})],$
(1)

where,

$\displaystyle
F(\tau_{\nu},\omega_{\nu})=\frac{1}{\exp(-\sqrt{3}\epsilon_{\nu}\tau_{\nu})(\epsilon_{\nu}-1)-(\epsilon_{\nu}+1)}\textsc{x}$
$\displaystyle\left[\frac{1-\exp(-(\sqrt{3}\epsilon_{\nu}+1/\mu)\tau_{\nu})}{\sqrt{3}\epsilon_{\nu}\mu+1}+\frac{\exp(-\tau_{\nu}/\mu)-\exp(-\sqrt{3}\epsilon_{\nu}\tau_{\nu})}{\sqrt{3}\epsilon_{\nu}\mu-1}\right],$
(2)

where $\omega_{\nu}=\frac{\sigma_{\nu}}{\kappa_{\nu}+\sigma_{\nu}}$ is the
albedo, which is defined by the scattering coefficient ($\sigma_{\nu}$) and
the absorption coefficient ($\kappa_{\nu}$), $\tau_{\nu}$ =
$\Sigma_{dust}\chi_{\nu}$, where $\chi_{\nu}={\kappa_{\nu}+\sigma_{\nu}}$,
$\mu=cos(i)$ considers the inclination effects by correcting the optical depth
by this inclination ($i$), and $\epsilon_{\nu}=\sqrt{1-\omega_{\nu}}$. This
model assumes that most of the dust content is settled in the mid-plane (see
Sierra et al., 2019 and Carrasco-González et al., 2019 for details). While
this assumption holds well for the outer HL Tau disk, it is not met within the
initial $\sim$ 32 au. A particle size distribution that follows a power-law
distribution was also assumed, specifically, $n(a)\propto a^{-p}$. We set the
value of the exponent $p$ to 3.5, which is a common value used in the
Interstellar Medium (ISM) (Mathis et al., 1977).

The dust opacity values were obtained assuming the DSHARP dust composition and
following the effective medium approximation (Birnstiel et al., 2018). As
usual in the study of protoplanetary disks, we first consider compact solid
spherical particles. However, we also perform an analysis by considering
spherical porous dust particles with a porosity of 90%. This is motivated by
very recent results which strongly suggest that the HL Tau disk contains
porous dust particles (Zhang et al. 2023 and Lin et al. 2023a).

The opacity tables for the porous particles were obtained using the same
procedure as (Birnstiel et al., 2018): the grain is assumed to be composed of
small monomers within a void matrix. The optical properties of the monomers
are computed first using the Bruggeman rule. This same rule is also used for
the compact grains. The monomers are then mixed with a void matrix using the
Maxwell-Garnett rule. It is well-established that the fluxes and,
consequently, dust masses exhibit significant dependence on the assumed dust
composition, particularly the quantity of carbon. For instance, when using the
opacities presented by Ricci et al. 2010 and Stadler et al. 2022 the fluxes
they inferred were higher which made the spectral indices of typical disks
more easily explained. Therefore, the composition of the grains plays a
crucial role, and the total dust mass derived in the subsequent section should
be approached with caution.

We model the SED of the disk at each radius independently and we fitted the
dust parameters (Tdust, amax and $\Sigma_{dust}$). We assume a flux
calibration uncertainty of 10% at 0.45 mm, 0.9 mm, 1.3 mm and 7 mm, and 5% at
2.1 mm. We compute the posterior probability distribution of our model at each
radius using the Markov Chain Monte Carlo (MCMC) implementation within emcee
(Foreman-Mackey et al., 2013). We assume uniform priors in dust surface
density and maximum grain size. Following Macías et al. 2021, we use a
conservative (i.e., wide) prior for the dust temperature based on the expected
radial temperature profile for a passively irradiated disk,
$T_{d}=(\phi{L}_{\star}/8\pi r^{2}\sigma_{\rm SB})^{0.25}$, assuming
${L}_{star}=6\pm 5\leavevmode\nobreak\ {L}_{\odot}$ (i.e., a very conservative
uncertainty on the stellar luminosity) and a range in $\phi$ (i.e., flaring
angle) between 0.005 and 0.3. This prior is wide enough to have no effect on
the posterior of the temperature at most radii, only helping avoid reaching
unphysically high dust temperatures at a few radii within the wide gap between
60 and 90 au, where the emission is optically thin and, therefore, the dust
temperature and dust surface density are highly degenerated.

In Figures 3 and 4, we show the dust parameters for the HL Tau disk obtained
from the fitting with the two different particles, solid and porous,
respectively. In Appendix A, we show the model radial intensity profiles at
each wavelength and their uncertainties, obtained from the median, 16th, and
86th percentiles of 500 random SEDs computed from the MCMC chains. As
anticipated, the dust temperature is very well constrained for both cases,
compact and porous particles. We also note that the surface density is well-
constrained in both cases. Moreover, our findings reveal that dust particles
in the HL Tau disk are relatively large, ranging from millimeters to
centimeters in size, these resemble well with models of dust evolution with
traps located at different distances (Pinilla et al., 2015).

Our models find a sudden, and likely unphysical, increase in particle size at
the 13 au gap. However, figures 7 and 8 show that our models are in fact
incapable of reproducing the observed emission at these radii. We believe the
reason for this is our limited angular resolution, paired with the fact that
the 13 au gap is highly depleted of large dust particles. As such, the SED at
high frequencies is dominated by the small grains present in the gap, while at
longer wavelengths the emission comes mostly from the beam smearing of the
emission at the 20 au ring. Higher angular resolution observations at long
wavelengths, only feasible with the upcoming ngVLA, are therefore necessary to
accurately constrain the dust properties at these radii.

A distinct gap at 70 au, wide and with almost no large particles, is evident
in 3 and 4, consistent with prior observations by Carrasco-González et al.
2019. Interestingly, an increase in dust surface density is then present at
the 84 au and 96 au rings, with grain sizes still remaining at a few hundreds
of microns. We note that an infalling streamer has been found to be
interacting with the HL Tau disk at these radii (Garufi et al., 2022). We
speculate that the small dust particles of this late infalling material might
be responsible for the creation of this dust ring, which would explain why the
particles at these radii have not yet grown up to mm/cm sizes. More detailed
modeling of the interaction between the disk and infalling streamers would be
needed to confirm this scenario.

Furthermore, clear distinctions emerge between the compact and porous particle
solutions. In general, the model using compact particles exhibits lower
temperatures and dust surface densities across the entire range of radii. The
particle sizes predicted by the porous dust composition appear to be slightly
smaller, although we note that they have a significantly wide uncertainty
range.

Given that our new 0.45 mm data have allowed us to better constrain the radial
dust temperature profile, we can more robustly compare the position of the
rings and gaps in HL Tau with the expected positions of the snowlines of some
of the most relevant volatiles in protoplanetary disks. We plot in Figures 3
and 4 the expected radii of the snowlines of H2O, NH3, CO2, H2S, and CO.
Except for the H2O snowline, the uncertainties in the sublimation temperatures
of these volatiles now dominate the uncertainty in their positions. These
large uncertainties still make it hard to extract any robust conclusion,
despite our improved temperature profile. Given the abundance of substructures
in HL Tau, these snowlines will always fall close to a ring or gap. This makes
it hard to confidently reject that snowlines are at the origin of at least a
few of the ring substructures in HL Tau, as previously proposed (Zhang et al.,
2015; Carrasco-González et al., 2019). However, many other substructures
appear to be far from the expected snowlines as in other protoplanetary disks
(Bae et al., 2023). This applies mostly to the outer substructures, between 65
and 98 au, but also some between 36 to 45 au. We therefore conclude that
mechanisms other than snowlines (e.g., MHD effects or planet-disk
interactions) must be responsible for most, if not all, the ring substructures
in HL Tau.

We infer a total dust masses of 2.1${}^{+1.88}_{-0.84}$ and
6.3${}^{+1.57}_{-1.04}$ MJ from the compact and porous particles,
respectively. The dust masses assuming porous particles are thus larger by a
factor of $\sim 3$, as previously reported by Zhang et al. 2023. Our
measurement of the dust mass with compact grains is two times higher than the
mass obtained in Carrasco-González et al. 2019 assuming compact particles with
the same composition as ours (1.04 MJ). The primary difference between their
analysis and ours is the inclusion of Band 9 in our sample. These additional
data have allowed us to more robustly constrain the dust temperature, which in
turn affects the rest of the modeling. We obtained lower dust temperatures
across all radii. The surface density is comparable at the inner radius but it
has a more rapid decline compared to their study. Additionally, our analysis
reveals more pronounced substructures in the surface density. Notably, we
observe larger dust particle sizes by approximately one order of magnitude
within 60 AU, but beyond this point, the particle sizes decrease more
significantly compared to the findings in Carrasco-González et al. 2019. These
larger particles inside 60 au have lower absorption opacities, which explains
why we obtain a higher dust mass than Carrasco-González et al. 2019.
Additionally, Booth & Ilee 2020 estimated a lower limit for the total gas mass
of the disk to be 209.5 MJ. This reconciles well with our results of the
compact dust particles if we consider a gas-to-dust ratio of a 100. However,
the elevated mass associated with porous dust particles from our analysis
might imply a reduced gas-to-dust ratio $\sim$ 33 when compared to the total
disk mass of 209.5 MJ. This suggests that the total gas mass of the disk might
indeed still be larger than 209.5 MJ.

We note that our assumption of a 1-D slab model might not be accurate in the
inner radius ¡ 32 AU, where the emission at 0.45 mm is so optically thick that
it could be tracing layers of the disk above the midplane (see section 4.2).
If all wavelengths are effectively probing adjacent and close layers in the
disk, indicating minimal variations in temperature and dust size distribution
between them, then the assumption could still be valid. The midplane of
protoplanetary disks is generally vertically isothermal, so it is most likely
that, if our observations are tracing significantly different layers of the
disk, this has only an effect on the maximum grain sizes traced at each
wavelength. Sierra & Lizano (2020) explored these effects in more detail and
found that the high optical depths in the innermost regions of disks could
result in incorrectly measuring the maximum grain size to be of a few hundreds
of microns. While we cannot discard that our modeling is affected by these
effects, the fact that we find maximum grain sizes of $\sim 5$ mm in the inner
region, instead of the 200-300 microns found by Sierra & Lizano (2020)
suggests that our modeling is not severely affected by the high optical
depths.

###  4.2 Azimuthally asymmetric emission in the inner disk

As mentioned in Section 3, the first ring, centered at 21 au, shows an
azimuthal asymmetry at 0.45 mm that becomes fainter at longer wavelengths and
we believe this is contingent upon the resolution and sensitivity of the
observations. This asymmetry appears as an increase in brightness at the
northeastern side of the ring, covering almost 180 degrees and centered on the
minor axis of the disk (see Fig. 1). The position of this asymmetry strongly
suggests that it is caused by the combination of the disk’s inclination and
its high optical depth at the observed wavelength (see Fig. 5). Specifically,
in the northeastern (NE) part of the disk, we are directly observing the
emission from the internal wall of the first ring, which is directly
irradiated by the central star, increasing its temperature. On the other hand,
in the southwestern (SW) part of the disk there is no direct line of sight
toward the wall, so the wall emission is blocked by the optically thick
material in the line of sight (see Figure 5).

Given that the HL Tau outflow is blue-shifted in the NE and red-shifted in the
SW (Lumbreras & Zapata, 2014), we can infer that the NE corresponds to the far
side of the disk, while the SW corresponds to the close side of the disk to
the observed, this explains the position of the asymmetry in that side of the
minor axis of the disk.

The asymmetry, in this case, is fundamentally similar to those found in Class
0/I objects (e.g Guerra-Alvarado et al. 2024,Lin et al. 2023b,Ohashi et al.
2023) where the dust particles did not have time to settle or are affected by
significant turbulence preventing the grains to fall into the midplane (larger
dust scale height), remaining optically thick even at longer wavelengths with
ALMA. On the other hand, this asymmetry is fundamentally different from the
large-scale asymmetries observed in other disks, where the effects become more
pronounced at longer wavelengths (e.g., IRS48 (van der Marel et al., 2015) and
HD142527 (Casassus et al., 2015) which is contrary to the observed behavior in
HL Tau. These asymmetries are commonly associated with dust traps where dust
can grow up to cm sizes.

Figure 5: Diagram explaining the scenario proposed for the origin of the
asymmetry seen in HL Tau at 0.45 mm.

In order to confirm that the asymmetry is caused by a geometric effect, we use
a toy model to reproduce the observed asymmetry. We use RADMC-3D (Dullemond et
al., 2012) to perform radiative transfer models at 0.45 mm while varying
parameters such as the disk dust mass and scale height. We then compare the
model emission with the observed data to gain insights into the disk’s
physical properties and vertical structure. In this specific radiative
modeling, a generic protoplanetary disk model was utilized within RADMC-3D.
The disk dust masses and scale height were adjusted to explore different
scenarios and to reproduce the observed asymmetry. The optool package (Dominik
et al., 2021) was employed to compute the DSHARP dust particle opacities using
amin as 0.050 $\mu$m and amax as 3mm. For this toy model, we are only
interested in the inner regions of the disk, so we consider only emission up
to a radius of 32 au (the location of the second gap). We assumed that all the
grains are at the same height and we set a small inner radius $R_{in}$=1 au. A
gap between 8 au and 17 au was included in the model to capture the effects of
the asymmetry at the outer edge of the gap.

The generic protoplanetary disk model uses a density distribution that goes
as:

$\displaystyle\rho(r,z)=\frac{\Sigma(r)}{H_{p}\sqrt{2\pi}}\exp(-\frac{z^{2}}{2H^{2}_{p}}),$
(3)

where $r$ is the distance to the star from the disk, $H_{p}$(r) is the dust
scale height of the disk and $\Sigma(r)$, is the dust surface density defined
as:

$\displaystyle\Sigma(r)=\Sigma_{0}(\frac{r}{r_{out}})^{-1},$ (4)

The dust scale height in the radiative modeling is defined as a power-law-like
dependence given by:

$\displaystyle H_{p}(r)=H_{100}(\frac{r}{100AU})^{1+\Psi},$ (5)

where $\Psi$ is the value of the flaring index set by default to 0.14 and
$H_{100}$, is the value of the scale height at a distance of 100 AU from the
central star.

Finally, a total stellar mass of 1.7 M${\odot}$ (Pinte et al., 2016), a
stellar luminosity of 6 $L_{\odot}$ and an effective temperature of 4395 $K$
(White & Hillenbrand, 2004) were assumed. In addition, the inclination angle
for HL Tau was set to 46.72∘ and the distance to 147 pc.

Several models were run by varying the dust scale height and the dust disk
mass within 32 au (M${}_{Disk}^{<32au}$) based on previous studies (Pinte et
al. 2016; Liu et al. 2017). We showed this parameters and their combination in
Table LABEL:table:RADMC-3D_parameters.

Table 1: Model parameters explored M${}_{Disk}^{<32au}$ | 1 MJ | 3 MJ | 5 MJ
---|---|---|---
H/R | 0.01 | 0.01 | 0.01
| 0.03 | 0.03 | 0.03
| 0.05 | 0.05 | 0.05
| 0.08 | 0.08 | 0.08
| 0.1 | 0.1 | 0.1

The model images were then convolved using a Gaussian beam with the same beam
as the 0.45 mm image. The final convolved model images that seem to better
reproduce the emission of the observation are shown in Figure 6. We also
extracted the deprojected emission from each pixel in a ring ranging from
0.12$\arcsec$ to 0.15$\arcsec$ and around 360∘ to capture the increased
emission in the asymmetry. We show the normalized azimuthal profiles from the
observation and all models in Figure 6. We note that our goal with this toy
models is to show that this type of asymmetries can be obtained at 0.45 mm and
that one could constrain the dust mass and scale height from these images. As
such, we do not try to reproduce the exact flux of our observations, since
that would require a much more detailed modeling varying more disk parameters.
The azimuthal profiles for 1, 3, and 5 Jupiter masses (Inner disk mass + ring
mass) are also shown in Figure 6.

Figure 6: RADMC-3D models of the inner HL Tau disk compared with the Band 9
observation. Top Figure: Convolved images of the 9 different models with
RADMC-3D that better reproduce the observation, corresponding to masses of 1
MJ, 3 MJ, and 5 MJ, paired with scale height values of 0.08, 0.05, and 0.03
$H_{p}$. The disk mass is defined as the dust mass within 32 au. Bottom
panels: Normalized azimuthal profiles considering the emission between 0.12
and 0.15 arcseconds in the HL Tau Band 9 observation and all the RADMC-3D
models. The lowermost panels depict the difference, measured in sigmas,
between the observation and the respective models

Our modeling shows that an asymmetry like the observed one in HL Tau can be
produced by the combined effects of disk orientation, high optical depth, and
larger scale heights at 0.45 mm. From Figure 6, it is evident that the dust
mass plays a significant role in reproducing the observed asymmetry at 0.45 mm
and in models with not a large enough scale height, the emission doesn’t
produce an asymmetry.

If the asymmetry arises from the scenario we outlined before, we find that a
disk dust mass (within 32 au) of approximately 1 $M_{J}$, combined with a
scale height H/R of 0.08, a disk dust mass of 3 $M_{J}$ with a scale height
H/R of 0.05, or a disk dust mass of 5 $M_{J}$ with a scale height H/R slightly
smaller than 0.05 best reproduces the observed asymmetry in the observations.
In other words, the higher the dust mass, the less dust scale height one needs
to reproduce an asymmetry. Because of this, obtaining an independent
estimation of the dust mass within 32 au would allow for a better constraint
on the scale height of the disk. We employ the mass obtained from our multi-
wavelength analysis in section 4.1, acknowledging its known limitations, to
give an initial approximation for the scale height in this specific region of
the disk. From this analysis, we obtained a disk dust mass inside the first 32
au of 0.86${}^{+1}_{-0.368}$ MJ. Therefore, if our proposed scenario is the
reason of the asymmetry, our derived mass would imply a scale height $\sim$¿=
0.08 H/R in that specific region of the disk. As seen in Figure 6, to more
effectively determine the scale height in subsequent studies we need a more
restrictive constraint on the disk dust mass together with a model that
reproduces the flux and the actual structure of the HL Tau disk.

We note that the scale height necessary to reproduce the observed asymmetry is
significantly larger than the estimates from radiative transfer modeling of
the 0.9-3 mm emission, H/R ¡ 0.01 (Pinte et al., 2016). One possible
explanation for this inconsistency is that, in Pinte et al. 2016, the focus
was primarily on gaps and rings in the outer disk, rather than those within 30
au as constrained in this paper. This result then suggests that the HL Tau
disk may exhibit greater vertical extent (thickness) in the inner regions
while maintaining a flatter profile in the outer regions. Indeed, radial
variations in scale height have been identified for the two rings of HD163296,
with the 67au ring being very thick (¿4au) and the 100au ring being very flat
(¡1au) (Doi & Kataoka 2021, Liu et al. 2022).

While there are limited constraints on the radial variation of scale height,
this outcome implies that turbulence in disks might be higher in the inner
regions compared to the outer regions. This aligns well with the high
accretion rates on some disks that cannot be explained by their current low
turbulence values. In this context, the magnetically driven accretion disk
models suggested by Delage et al. 2022 could explain the heightened turbulence
parameter and increased scale height, leading to the observed asymmetry beyond
the dead zone (at 23 AU) and within the MRI-active layer in the HL Tau disk.
Additionally, this is consistent with what is expected in the case of vertical
shear instability with dust coagulation, where the transition between thick
and thin occurs at around 100 au (Pfeil et al. 2023). Finally, it also remains
inconclusive whether the streamer observed by Garufi et al. (2022) in HL Tau
could have an effect at smaller radii through the flow of material into the
inner disk. The findings in Jiang et al. (2024) suggest a potential influence
of this streamer on the turbulence at the outer radius (70 au), raising
questions about a possible impact on smaller scales. Further analysis is
essential to discern the factors contributing to the increased turbulence, if
any, in the HL Tau inner disk and unravel the specifics of this phenomenon.

###  4.3 Substructure in the inner disk: a sign of the traffic-jam effect?

As mentioned, our highest angular resolution image at 0.45 mm revealed an
increase in the brightness temperature within $\sim$2.5 au (see the inset in
Figure 2). There are several possible explanations for this.

First, we discard that this increase is due to free-free contamination. It is
well known that HL Tau drives a powerful jet whose emission has been
previously modeled using high angular resolution VLA data at long wavelengths
(Carrasco-González et al., 2019). Thus, we know that in the 23-43 GHz
frequency range, the expected free-free emission from the radio jet is given
by

$\displaystyle\left[\frac{F_{\nu}}{\mu Jy\
{beam}^{-1}}\right]=175\times\left[\frac{\nu}{32.5\ GHz}\right]^{0.7}$ (6)

where $\nu$ is the frequency of the image at which we want to calculate the
free-free emission. Note that in this frequency range, the free-free emission
was found to be partially optically thin, which is reflected in the 0.7
spectral index. Now, since the optical depth of the free-free emission is
proportional to $\nu^{-2.1}$ (e.g., Reynolds 1986), at the frequency of our
0.45 mm image, i.e. $\sim$650 GHz, we expect free-free emission from the jet
to be optically thin. Thus, there must be a turning point between 50 and 650
GHz. A conservative assumption is that emission from the jet becomes totally
optically thin at 100 GHz. Then, using equation 6 to estimate the free-free
emission at 100 GHz, and extrapolating this with a spectral index of $-$0.1 to
higher frequencies, we obtain that an estimate of the free-free contamination
at 650 GHz would be $\sim$300 $\mu$Jy, which is a $\sim$2% of the peak
emission at the center of the disk in the 0.45 mm image. Note that the actual
contribution at this frequency should be lower since it is very likely that
the free-free emission would become optically thin at a frequency lower than
100 GHz. Another potential explanation for a significant free-free emission at
the central part of the disk is the possibility of it becoming highly
optically thick due to photoionization. However, this scenario faces
challenges because HL Tau is a K-type star, and such stars typically do not
emit sufficient ionizing photons to fully ionize a very dense gas. Taking all
of these factors into account, it leads us to the conclusion that the elevated
emission in the inner region of the disk is unlikely to be primarily caused by
free-free contributions.

An interesting explanation for the increased dust emission at the inner part
of the disk is the possibility of a dust particle pile-up occurring in that
region. In fact, our unprecedented angular resolution allows us to probe the
dust emissions down a few au, inside the estimated water snowline at $\sim$ 5
au if one assumes $T_{disk}=T_{B}$ (Zhang et al., 2015, Figure 2). Inside the
water snowline, theoretical models predict the appearance of a traffic-jam
effect (e.g Banzatti et al. 2015, Pinilla et al. 2017), as dry silicates are
more sensitive to fragmentation than ice-rich grains (Dominik & Tielens 1997;
Supulver et al. 1997; Wada et al. 2013; Gundlach & Blum 2014), leading to
lower particle size and reduced drift speed. The traffic-jam effect offers a
favored pathway for the formation of dry planetesimals in the inner disk,
however, it has been recently put into question in the light of new
experiments on the stickiness of ice-rich grains (e.g. Gundlach et al. 2018;
Musiolik & Wurm 2019) and observations of the resilience of icy pebbles to
sublimation (Houge et al., 2023). Interestingly, Facchini et al. (2024)
recently detected water emission at 321 GHz coming from the inner disk of HL
Tau, but their spatial resolution was not enough to robustly measure the
position of the water snowline.

Our results could be consistent with the presence of a traffic-jam effect, as
such a pile-up of solids would increase the dust density and optical depth
inside the water snowline (e.g. Banzatti et al. 2015), causing observations to
probe higher layers of the disk with more elevated dust temperatures, thereby
manifesting a higher brightness temperature as we see in Figure 2 (See
Drażkowska & Alibert 2017 and Ros et al. 2019). Moreover, an increase in the
population of smaller particles will also result in a decrease in the albedo,
which in turn results in an increase in the emission.

A second possibility is that the increase in the brightness temperature is
attributed to viscous heating occurring in the inner disk. Overall, the
heating mechanisms in protoplanetary disks significantly influence the
temperature profile of these. The two dominant heating processes in these
disks are stellar radiation and viscous heating (D’Alessio et al., 2005).
Across the majority of the disk’s extent, stellar radiation contributes the
most to the disk’s thermal structure. However, in the innermost regions of the
disk close to the mid-plane, where dust surface densities reach notably high
values, it is highly likely that viscous heating becomes the primary
contributor to the disk temperature profile (D’Alessio et al., 2001).More
recently, Chambers (2009) found that the mid-plane temperature may be affected
by viscous heating within the radial range of 1 to 40 au in protoplanetary
disks. This phenomenon also influences the position of the snowline, which
varies depending on the age of the disk. The position of these, however, is
dependent upon the characteristics of the star and many other parameters which
are hard to constrain.

As of now, it is challenging to distinguish between the traffic-jam or viscous
heating scenarios, but we strongly encourage future works to focus on that
aspect given the implications it could have concerning planet formation in the
inner disk. For example, high angular resolution observations at other
wavelengths could help constrain the dust emission spectral index and particle
size, which would help in understanding the underlying mechanisms driving the
observed changes in brightness temperature.

##  5 Summary and conclusions

We report new sensitive ALMA images at 0.45 mm of the HL Tau disk with an
unprecedented high angular resolution of 12 mas or 1.7 au. At this wavelength,
the detected emission is known to be highly optically thick. This is evidenced
by the shallower substructures observed in our observation. The extended
nature of the disk suggests that we are tracing smaller dust particles that
are more closely coupled to the gas.

By combining the Band 9 image with the images at other wavelengths, we have
updated the model of the dust properties in the disk. The 0.45 mm image has
allowed us to better constrain the dust temperature of the disk, which
provided better constraints on the dust density and particle size. The
particle sizes are still predicted to be larger than mm sizes until $\sim 60$
au but there is a significant decrease (to $\sim$ 200 $\mu$m) between 60-100
au. These smaller particles coincide with the region where an accretion shock
of infalling material has been found traced by SO and SO2. The final derived
disk dust masses are 2.1 MJ for compact particles and 6.3 MJ for porous dust
particles. These dust masses translate into gas-to-dust ratios $\leq$ 100.

Our 0.45 mm data clearly shows an asymmetry in the emission coming from the
first ring. The emission is considerably more intense in the NE part of the
disk than in the SW part. We interpreted this as a combination of very
optically thick emission, a moderately inclined disk, and a large dust scale
height. In the NE part of the disk, we are directly viewing the emission from
the internal wall of the first ring, which is illuminated by the radiation
from the central star. This interpretation is supported by our radiative
transfer model which demonstrates the dependence of this asymmetry with the
dust mass of the inner regions of the disk and the dust scale height. We
propose that this kind of asymmetries in other disks can be used to constrain
the dust scale height if an independent estimation of the dust mass of the
disk has been made. Based on our proposed scenario and the multi-wavelength
analysis model, we implied the possibility of a dust scale height H/R $\sim$
¿= 0.08 for the HL Tau disk at least in the inner 32 au. This approximation
seems in large contrast with the constraints identified by Pinte et al. 2016
who determined H/R ¡ 0.01 in the outer disk. We note that this discrepancy
might stem from higher turbulence in the inner regions of HL Tau compared to
its outer region leading to radial variations of the vertical structure of the
HL Tau disk. Incorporating vertical structure in the multi-wavelength analysis
to determine the disk mass and employing more quantitative modeling
approaches, rather than qualitative ones in RADMC-3D, may be necessary to
obtain a more realistic scale height of the disk.

We have resolved the innermost part of the HL Tau disk showing that there is
some previously unseen substructure, revealed as a steep increase of the
brightness temperature very close to the central star, at radii smaller than 5
au. Interestingly, the radius at which this increase starts appears to
coincide with the position of the water ice line. Therefore, we discuss the
possibility of a pile-up of small particles due to a decrease in the drift
velocity of the small grains after crossing the water snowline and suffer a
more efficient fragmentation. We also discuss a different possibility in which
we are resolving the region where the heating of the disk changes from being
dominated by irradiation to being dominated by viscous heating. Previous
theoretical models predict that this change should take place very close to
the protostar. To further confirm either of these scenarios, additional
modeling and observations at other wavelengths are necessary for the HL Tau
disk.

###### Acknowledgements.

This paper makes use of the following ALMA data:
ADS/JAO.ALMA$\\#$2011.0.00005.E and ADS/JAO.ALMA$\\#$2017.1.01178.S. ALMA is a
partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. C.C.-G. acknowledges
support from UNAM DGAPA-PAPIIT grant IG101321 and from CONACyT Ciencia de
Frontera project ID 86372. The research of M.V. was supported by an
appointment to the NASA Postdoctoral Program at the NASA Jet Propulsion
Laboratory, administered by Oak Ridge Associated Universities, under contract
with NASA.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * Asaki et al. (2020) Asaki, Y., Maud, L. T., Fomalont, E. B., et al. 2020, ApJS, 247, 23
  * Bae et al. (2023) Bae, J., Isella, A., Zhu, Z., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 423
  * Banzatti et al. (2015) Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15
  * Barrière-Fouchet et al. (2005) Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45
  * Booth & Ilee (2020) Booth, A. S. & Ilee, J. D. 2020, MNRAS, 493, L108
  * Carrasco-González et al. (2016) Carrasco-González, C., Henning, T., Chandler, C. J., et al. 2016, ApJ, 821, L16
  * Carrasco-González et al. (2019) Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71
  * Casassus et al. (2015) Casassus, S., Wright, C. M., Marino, S., et al. 2015, ApJ, 812, 126
  * Cazzoletti et al. (2018) Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161
  * Chambers (2009) Chambers, J. E. 2009, ApJ, 705, 1206
  * D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321
  * D’Alessio et al. (2005) D’Alessio, P., Calvet, N., & Woolum, D. S. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 341, Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth, 353
  * de Juan Ovelar et al. (2013) de Juan Ovelar, M., Min, M., Dominik, C., et al. 2013, A&A, 560, A111
  * Delage et al. (2022) Delage, T. N., Okuzumi, S., Flock, M., Pinilla, P., & Dzyurkevich, N. 2022, A&A, 658, A97
  * Dipierro et al. (2016) Dipierro, G., Laibe, G., Price, D. J., & Lodato, G. 2016, MNRAS, 459, L1
  * Doi & Kataoka (2021) Doi, K. & Kataoka, A. 2021, ApJ, 912, 164
  * Dominik et al. (2021) Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library, ascl:2104.010
  * Dominik & Tielens (1997) Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647
  * Drażkowska & Alibert (2017) Drażkowska, J. & Alibert, Y. 2017, A&A, 608, A92
  * Dullemond & Dominik (2004) Dullemond, C. P. & Dominik, C. 2004, A&A, 421, 1075
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015
  * Facchini et al. (2024) Facchini, S., Testi, L., Humphreys, E., et al. 2024, Nature Astronomy
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Furlan et al. (2008) Furlan, E., McClure, M., Calvet, N., et al. 2008, ApJS, 176, 184
  * Galli et al. (2018) Galli, P. A. B., Loinard, L., Ortiz-Léon, G. N., et al. 2018, ApJ, 859, 33
  * Garufi et al. (2022) Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104
  * Guerra-Alvarado et al. (2024) Guerra-Alvarado, O. M., van der Marel, N., Di Francesco, J., et al. 2024, A&A, 681, A82
  * Guidi et al. (2022) Guidi, G., Isella, A., Testi, L., et al. 2022, A&A, 664, A137
  * Gundlach & Blum (2014) Gundlach, B. & Blum, J. 2014, in Asteroids, Comets, Meteors 2014, ed. K. Muinonen, A. Penttilä, M. Granvik, A. Virkki, G. Fedorets, O. Wilkman, & T. Kohout, 191
  * Gundlach et al. (2018) Gundlach, B., Schmidt, K. P., Kreuzig, C., et al. 2018, MNRAS, 479, 1273
  * Hashimoto et al. (2023) Hashimoto, J., Liu, H. B., Dong, R., et al. 2023, AJ, 166, 186
  * Houge et al. (2023) Houge, A., Macías, E., & Krijt, S. 2023, arXiv e-prints, arXiv:2312.01856
  * Huang et al. (2018) Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJ, 869, L43
  * Isella et al. (2016) Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101
  * Jiang et al. (2024) Jiang, H., Macías, E., Guerra-Alvarado, O. M., & Carrasco-González, C. 2024, A&A, 682, A32
  * Jin et al. (2016) Jin, S., Li, S., Isella, A., Li, H., & Ji, J. 2016, ApJ, 818, 76
  * Kataoka et al. (2017) Kataoka, A., Tsukagoshi, T., Pohl, A., et al. 2017, ApJ, 844, L5
  * Lin et al. (2023a) Lin, Z.-Y. D., Li, Z.-Y., Stephens, I. W., et al. 2023a, arXiv e-prints, arXiv:2309.10055
  * Lin et al. (2023b) Lin, Z.-Y. D., Li, Z.-Y., Tobin, J. J., et al. 2023b, ApJ, 951, 9
  * Liu et al. (2022) Liu, Y., Bertrang, G. H. M., Flock, M., et al. 2022, Science China Physics, Mechanics, and Astronomy, 65, 129511
  * Liu et al. (2017) Liu, Y., Henning, T., Carrasco-González, C., et al. 2017, A&A, 607, A74
  * Long et al. (2020) Long, F., Pinilla, P., Herczeg, G. J., et al. 2020, ApJ, 898, 36
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
  * Lumbreras & Zapata (2014) Lumbreras, A. M. & Zapata, L. A. 2014, AJ, 147, 72
  * Macías et al. (2021) Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Musiolik & Wurm (2019) Musiolik, G. & Wurm, G. 2019, ApJ, 873, 58
  * Ohashi et al. (2023) Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8
  * Ohashi et al. (2020) Ohashi, S., Kataoka, A., van der Marel, N., et al. 2020, ApJ, 900, 81
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519
  * Pfeil et al. (2023) Pfeil, T., Birnstiel, T., & Klahr, H. 2023, arXiv e-prints, arXiv:2310.07332
  * Pinilla et al. (2012a) Pinilla, P., Benisty, M., & Birnstiel, T. 2012a, A&A, 545, A81
  * Pinilla et al. (2012b) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012b, A&A, 538, A114
  * Pinilla et al. (2015) Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015, A&A, 573, A9
  * Pinilla et al. (2021) Pinilla, P., Lenz, C. T., & Stammler, S. M. 2021, A&A, 645, A70
  * Pinilla et al. (2017) Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
  * Pizzati et al. (2023) Pizzati, E., Rosotti, G. P., & Tabone, B. 2023, MNRAS, 524, 3184
  * Rau & Cornwell (2011) Rau, U. & Cornwell, T. J. 2011, A&A, 532, A71
  * Reynolds (1986) Reynolds, S. P. 1986, ApJ, 304, 713
  * Ricci et al. (2010) Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15
  * Ros et al. (2019) Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65
  * Sierra & Lizano (2020) Sierra, A. & Lizano, S. 2020, ApJ, 892, 136
  * Sierra et al. (2019) Sierra, A., Lizano, S., Macías, E., et al. 2019, ApJ, 876, 7
  * Sierra et al. (2021) Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14
  * Stadler et al. (2022) Stadler, J., Gárate, M., Pinilla, P., et al. 2022, A&A, 668, A104
  * Supulver et al. (1997) Supulver, K. D., Bridges, F. G., Tiscareno, S., Lievore, J., & Lin, D. N. C. 1997, Icarus, 129, 539
  * Tazzari et al. (2021) Tazzari, M., Testi, L., Natta, A., et al. 2021, MNRAS, 506, 5117
  * van der Marel et al. (2015) van der Marel, N., Pinilla, P., Tobin, J., et al. 2015, ApJ, 810, L7
  * van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199
  * Villenave et al. (2020) Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164
  * Villenave et al. (2022) Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11
  * Wada et al. (2013) Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62
  * White & Hillenbrand (2004) White, R. J. & Hillenbrand, L. A. 2004, ApJ, 616, 998
  * Zhang et al. (2015) Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7
  * Zhang et al. (2023) Zhang, S., Zhu, Z., Ueda, T., et al. 2023, ApJ, 953, 96

##  Appendix A MCMC model vs observations: VLA fitting

Figure 7: Comparison of the brightness temperatures of the model (dash lines)
with the observations (solid lines) across each radius and at each wavelength.
In this Figure, the compact dust particles were used, and the VLA emission is
plotted and fitted within the MCMC. The vertical lines are the positions of
the gaps and bright rings taken from the dust parameters and the size of the
beam is plotted as a Gaussian on the top right corner for the 5 fitted
wavelengths (50 mas). Figure 8: Same as Figure B.1 but with the porous
particle model instead.


# Geometrical analysis of Kerr-lens mode-locking for high-peak-power ultrafast
oscillators

Peiyu Xia1, Makoto Kuwata-Gonokami1, Kosuke Yoshioka${}^{1,2^{*}}$ 1Department
of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan  2Photon Science Center, School of
Engineering, The University of Tokyo, 2-11-16 Yayoi, Bunkyo-ku, Tokyo
113-8656, Japan<EMAIL_ADDRESS>
###### Abstract

Ultrashort pulses from Kerr-lens mode-locked oscillators have inspired a
variety of applications. The design and alignment of these laser resonators
have thus far been theoretically supported by the conventional analysis of
beam propagation. However, the well-established theoretical framework is
sometimes beyond the scope of high-peak-power oscillators. In this paper, we
analyze the geometry of ring resonators by extending the ABCD-matrix method to
a high-peak-power regime. The guidelines to achieving stable Kerr-lens mode-
locking is provided for high-peak-power pulses.

  * March 2020

##  1 Introduction

Following the demonstration of Kerr-lens mode-locking (KLM) by using the
optical Kerr effect in Ti:Sapphire [1], laser resonators with a Kerr medium
have become the basic platform for the generation of ultrashort optical
pulses. Besides, the emergence of Ti-, Yb-, and Cr-doped gain crystals allows
to generate femtosecond pulses in the visible to the near-infrared region [2,
3, 4, 5]. In designing these Kerr-lens mode-locked oscillators, one of the
essential tasks is to judge if the laser can establish stable KLM.

To evaluate the stable KLM operation quantitatively, earlier theoretical works
have presented a basic framework described as the following [6, 7]. First, KLM
is caused by the power-dependent gains or losses induced by the self-focusing
effect in a Kerr medium. Second, the slope of the gains or losses is evaluated
by numerically analyzing the beam propagation in the spatial domain. Here, the
essential parameter to evaluate the slope is the relative spot-size variation
defined as follows:

$\eqalign{\delta=\left[-\frac{1}{\omega}\frac{d\omega}{dP}\right]_{P=0}.}$ (1)

$\omega(P)$ is the spot size at the Kerr medium as a function of the
intracavity power, $P$. So far, the $\delta$ parameter has been mostly used as
a measure of stability to provide the basic guidelines for KLM in the
commonly-used linear or ring resonators, and to design these resonators [8, 9,
10, 11, 12, 13, 14]. Although the more realistic frameworks which treat the
spatiotemporal dynamics have been developed [15, 16], the analysis based on
the $\delta$ parameter is still widely used [17, 18, 19, 20, 21] due to its
simplicity and low cost of computation.

However, the $\delta$ parameter is not always valid when the intracavity power
$P$ is close to or larger than the critical power for self-focusing,
$P_{\mathrm{c}}$, defined as

$P_{\mathrm{c}}=\frac{\alpha\lambda^{2}}{4\pi n_{0}n_{2}},$ (2)

where $n_{0}$ and $n_{2}$ are the linear and nonlinear refractive indices of
the Kerr medium, $\lambda$ is a laser wavelength, and $\alpha$ is a constant
calculated to be 1.86–2.0 for Gaussian beam [7, 22]. The theoretical work by
Larotonda showed that the stability of KLM saturates when $P/P_{\mathrm{c}}$
approaches 1 in a Kerr medium [23]. Besides, Bartels and Kurz pointed out that
experimentally realized oscillators including their result usually satisfied
$P/P_{\mathrm{c}}>1$ [24]. For instance, the Ti:Sapphire resonator of Ref. [1]
satisfied $P/P_{\mathrm{c}}=1.3$, calculated from the reported output peak
power of 90 kW with a 3.5% output coupler and $P_{\mathrm{c}}=2$ MW for
Ti:Sapphire. In the same way, the Cr4+:YAG resonator of Ref. [3] satisfied
$P/P_{\mathrm{c}}=1.4$, and the Cr2+:ZnSe resonator of Ref. [5] satisfied
$P/P_{\mathrm{c}}=1.2$. We note that these works are the first demonstration
of KLM using the corresponding gain crystals. Therefore, it is practically
necessary to understand the mechanism of stable KLM in laser resonators with
$P/P_{\mathrm{c}}>1$.

In this work, we numerically analyze beam propagation in ring resonators
particularly in the high-peak-power regime ($P/P_{\mathrm{c}}$ up to 1.6).
Based on the ABCD-matrix formulation, we find that the stable region of the
ring resonator for the pulse propagation differed significantly from that for
the continuous-wave (CW) propagation which corresponds to
$P/P_{\mathrm{c}}=0$. Moreover, we extend the analysis based on the $\delta$
parameter to better evaluate the stability of KLM.

##  2 Method

Figure 1: (a) Configuration of a typical Kerr-lens mode-locked ring
oscillator, and (b) the representative model for the ABCD-matrix calculation.
The black dotted line represents the periodic propagation of a circulating
beam and corresponds to the fixed position for the calculation of the spot
size. $x$ is the displacement of the medium from the center position, and $d$
is the thickness of the medium, neglecting the oblique incidences, $z$ and
$L_{\mathrm{c}}$, are the shorter and longer propagation distances between the
concave mirrors, CM1 and 2, respectively. The inset illustrates the assumption
of the pump-to-laser overlap and the position to calculate $\delta_{K}$. L1
and 2: a pair of lenses with the same focal length as CM1 and 2.

As shown in Fig. 1, we analyze the ring resonator with neglect of the
astigmatism of oblique incidence to a Kerr medium and mirrors. Consequently,
the beam propagation in the ring resonator is equally represented as in Fig.
1(b). We assume that the spot size of the pump beam is always aligned smaller
than that of the resonant intracavity beam, as shown in the inset of Fig.
1(b). This alignment is a common technique called “soft-aperture KLM” [6],
where the higher gain favors the smaller spot size in the gain medium. The
radius of curvature of each concave mirror is 30 mm, corresponding to the
focal length, $f$, of 15 mm. The thickness, $d$, and refractive index,
$n_{0}$, of the medium is 2.5 mm and 1.76, respectively. $L_{\mathrm{c}}$ is
180 mm, and the laser wavelength, $\lambda$, is 800 nm. These are the
conventional parameters used for KLM in the Ti:Sapphire ring resonators with
GHz-level repetition rates [24, 25], which are also estimated to satisfy
$P/P_{\mathrm{c}}>1$. Compared with the linear resonator, the ring resonator
has a single stable region, as described in Eq. 3, whereas the linear
resonator has double stable regions [6]. Therefore, the understanding of KLM,
based on the ring resonator, is clearer and analogous to the linear resonator.

We first use the formalism based on the ABCD matrix to calculate the self-
consistent Gaussian mode. In the case of the CW propagation, the stable region
is given by

$2f+d-d_{\mathrm{e}}<z<\frac{2fL_{\mathrm{c}}}{L_{\mathrm{c}}-2f}+d-d_{\mathrm{e}},$
(3)

where $d_{\mathrm{e}}=d/n_{0}$ is the effective length of the medium, and $z$
is the distance between the concave mirrors. In our case, Eq. 3 corresponds to
$31.08\ \mathrm{mm}<z<37.08\ \mathrm{mm}$. These lower and upper bounds are
defined as $z_{\mathrm{min}}$ and $z_{\mathrm{max}}$ repectively and indicated
in Fig. 2. In the case of the pulse propagation, nonlinear propagation in the
Kerr medium caused by self-focusing is expressed by the following matrix [7].

$\mathbf{M}(\gamma)=\sqrt{1-\gamma}\left(\begin{array}[]{cc}1&d_{\mathrm{e}}\\\
-\frac{\gamma}{(1-\gamma)d_{\mathrm{e}}}&1\end{array}\right).$ (4)

The $\gamma$ parameter depends on the intracavity power and is expressed by

$\gamma=K\left[1+\frac{1}{4}\left(\frac{2\pi\omega_{\mathrm{c}}^{2}}{\lambda
d_{\mathrm{e}}}-\frac{\lambda
d_{\mathrm{e}}}{2\pi\omega_{0}^{2}}\right)^{2}\right]^{-1},$ (5)

where $K=P/P_{\mathrm{c}}$ is the normalized intracavity power,
$\omega_{\mathrm{c}}$ is the spot size at the center of the medium, and
$\omega_{0}$ is the spot size at the beam waist calculated as if $K=0$. The
formalism using both linear ABCD matrices and the matrix of Eq. 4 allows to
calculate the spot size for the pulse propagation. However, Eq. 5 has
difficulty with multivariable analysis, because both $\omega_{\mathrm{c}}$ and
$\omega_{0}$ depend on the resonator parameters and the intracavity peak
power. For multivariable analysis, we modify Eq. 5 as follows.

$\gamma=K\left[1+\left(\frac{d^{-1}+\mathrm{Re}\left(q^{-1}\right)}{\mathrm{Im}\left(q^{-1}\right)}\right)^{2}\right]^{-1}.$
(6)

where $q^{-1}=R_{\mathrm{s}}^{-1}-i\lambda/\pi n_{0}\omega_{\mathrm{s}}^{2}$
is the reciprocal of $q$, given by the radius of the curvature,
$R_{\mathrm{s}}$, and the spot size, $\omega_{\mathrm{s}}$, at the incident
surface of the medium. After the parameters are set, the calculation of the
round-trip ABCD matrix of the resonator, including matrix
$\mathbf{M}(\gamma)$, leads us to find the $q$ parameter of any plane, which
becomes the iterative solution to the self-consistent equation,
$q=[A(\gamma)q+B(\gamma)]/[C(\gamma)q+D(\gamma)]$, with $\gamma=\gamma(q)$
given by Eq. 6.

Next, we define the following $\delta_{K}$ parameter, extended from Eq. 1, to
evaluate the stability of KLM.

$\delta_{K}=\left[-\frac{1}{\omega}\frac{d\omega}{d\left(P/P_{\mathrm{c}}\right)}\right]_{P=KP_{\mathrm{c}}}.$
(7)

Eq. 7 represents the relative spot size variation in perturbation to the
normalized peak power of $K$, i.e. the robustness against the spatial
perturbation. The case of $K=0$ corresponds to Eq. 1. The negatively larger
$\delta_{K}$ in the gain medium leads to a higher gain under the soft-aperture
condition as represented in the inset of Fig. 1(b). Compared with Eq. 1, the
$\delta_{K}$ parameter deviates from the $\delta$ parameter when $K$
approaches 1. Therefore, the evaluation of $\delta_{K}$ is more desirable to
achieve intracavity peak power around $P_{\mathrm{c}}$ during the KLM process.
For example, if $\delta_{K}$ is positive at certain peak power, the pulse
propagation is unstable due to no gain or losses, which results in a decrease
of the peak power. This gain-loss balance also causes the mode-locking process
in the time domain, as described in the role of saturable absorbers and self-
amplitude modulation [26, 27, 28].

##  3 Result

Figure 2: Variation of the spot size at the fixed position as indicated in
Fig. 1(b). The normalized intracavity peak power, $K$, is (a) 0.0, (b) 0.4,
(c) 0.8, and (d) 1.6, respectively. The purple arrows indicate the inner and
outer stability edges for the CW propagation expressed by Eq. 3.

Figure 2 shows the variation of the spot size as a function of $z$ and $x$,
the position of which is fixed in the resonator (see Fig. 1(b)). We found the
self-consistent solution only in the colored region, where the calculated spot
size is finite. The finite value out of the color range is also displayed with
the maximum or minimum color. In some conditions near the stability edge, we
found two self-consistent solutions, as reported on the linear resonator [7].
However, one is always discontinuous as a function of $z$ and $x$ and is
therefore excluded. Figure 2(a) shows that the spot size for the CW
propagation ($K=0$) is independent on the medium position $x$ at any $z$. On
the other hand, significant $x$-dependencies are observed in Figs. 2(b)–(d).
When the intracavity power is $0.4P_{\mathrm{c}}$ as in Fig. 2(b), the spot-
size variation changes significantly from that for the CW propagation,
especially near the stability edge ($z\simeq z_{\mathrm{min}}$ and $z\simeq
z_{\mathrm{max}}$). When the intracavity power is higher than $P_{\mathrm{c}}$
as is the case of Fig. 2(d), parts of the self-consistent solution are no
longer found near the inner stability edge with the medium close to the center
position ($z\simeq z_{\mathrm{min}}$ and $x\simeq 0$). This feature indicates
that the initially stable CW propagation can evolve into the unstable pulse
propagation during the KLM process. The contrasting feature in Figs. 2(b)-(d)
is the appearance of two pairs of sharp edges in the stable region with
increasing the intracavity power. The sharp edges grow into the unstable
region for the CW propagation ($z<z_{\mathrm{min}}$ or $z_{\mathrm{max}}<z$),
as can be seen from Fig. 2(a) to Fig. 2(d). This behaviour indicates that the
stable region for the pulse propagation can extend well beyond that for the CW
propagation, and finally into a fish-like shape. Therefore, the initially
unstable CW propagation can be developed into the stable pulse propagation,
provided that the concave mirrors and the medium position, $z$ and $x$, are
properly aligned. These differences of the stable region between the pulse and
CW propagation provide insight into some of the experimental results that
showed the increase of the laser power by an order of magnitude during the KLM
process, such as 1.2 W from 20–40 mW [24], and 900 mW from $\sim$50 mW [25].

Figure 3: Variation of the stability parameter, $\delta_{K}$. The normalized
intracavity peak power, K, is (a) 0.0, (b) 0.4, (c) 0.8, and (d) 1.6,
respectively. The contour lines enclose the areas having negatively large
values.

Next, we evaluate the pulse stability to find the optimal condition to achieve
stable KLM in the high-peak-power regime. Figure 3 shows the variation of
$\delta_{K}$ as a function of the same resonator parameters as used in Fig. 2.
In each of the resonator conditions, the $\delta_{K}$ parameters are
calculated at the position of the minimum spot size inside the gain medium, as
indicated in the inset of Fig. 1(b). When the intracavity power is much lower
than $P_{\mathrm{c}}$ ($K\ll 1$) which almost corresponds to Fig. 3(a), the
stability parameter is almost the same as given by Eq. 1. The negatively large
value can be found near the inner stability edge with the medium close to the
center position ($z\simeq z_{\mathrm{min}}$ and $x\simeq 0$). This condition
of alignment is consistent with the experimental condition needed to achieve
soft-aperture KLM of the ring resonators [30,31]. In contrast, the high-peak-
power pulses do not favor this condition, because the self-consistent solution
will disappear with increasing the intracavity power (see from Fig. 2(a) to
Fig. 2(d)). Instead of this unfavorable condition, if the medium is displaced
from the center position to some extent, the self-consistent solution can be
found near the sharp edges (ex. $1.5\ \mathrm{mm}\leq|x|$ for $K=1.6$ as in
Fig. 2(d)). Another item of note is the condition where $\delta_{K}$ has a
negatively large value as mapped in Fig. 3. The negatively large area (deep
blue area) is regarded as optimal to obtain the pulses with the respective
peak power, because the negatively large value corresponds to the amount of
the pulse stability. These areas are mainly found at $z\simeq
z_{\mathrm{min}}$ and $|x|\simeq 1.2,\ 1.5,\ \mathrm{and}\ 2.0$ mm for Figs.
3(b)–(d), respectively. Therefore, the proper displacement of the Kerr medium
according to the intracavity power is essential to achieve stable KLM.

We note that the other negatively large areas are also found in the unstable
region for the CW propagation ($z<z_{\mathrm{min}}$ or $z_{\mathrm{max}}<z$).
These conditions may allow pulses to be built up during the KLM process by a
strong mechanical perturbation to a laser resonator, or by an external
injection seeding.

Figure 4: (a) Power scaling in soft-aperture KLM with respect to the medium
displacement. The stability parameter $-\delta_{K}(x)$ is normalized to 1 for
each value of $K$. (b)-(d) The pulse propagation (red lines) around the medium
is simulated on the respective conditions as marked in (a). The CW propagation
(blue dashed lines) is shown at the same displacement assuming linear
propagation (K=0).

To understand the dynamics of KLM with respect to the medium position, figure
4 shows the three types of behaviour in the pulse propagation caused by the
strong self-focusing effect inside the medium. As shown in Fig. 4(a), we fix
the $z$ parameter to 31.2 mm, which is close to the stability edge
($z_{\mathrm{min}}$), and analyze the pulse propagation with the fixed power
of $2.0P_{\mathrm{c}}$. First, figure. 4(b) indicates that the beam with its
focus position inside the medium causes the self-focusing collapse. This
condition is satisfied almost when the displacement of the medium is shorter
than the effective thickness, $d_{\mathrm{e}}=d/n_{0}$, The collapse of the
beam could not result in stable KLM, due to damage to the medium or pulse
breakup. Second, when the displacement of the medium is near the optimal area
(deep blue area), the focus position of the beam is found to be usually close
to the surface of the medium. In this condition as shown in Fig. 4(c), the
counterbalance between linear diffraction and self-focusing results in the
effective self-trapping effect. This behaviour suggests that the self-trapping
effect contributes significantly to the stable KLM in the high-peak-power
regime, and moreover to the appearance of the sharp edges in Fig. 2(b)–(d). At
last, the further displacement of the medium result in an approach to the CW
mode as shown in Fig. 4(d). Although the pulse propagation may be stable due
to the relatively small but negative value of $\delta_{K}$, the condition of
Fig. 4(c) may as well be employed.

To summarize the discussion in Figs. 2–4, the following guidelines are
provided to achieve soft-aperture KLM in ring resonators.

  * •

The alignment condition in the CW operation should be close to the inner
stability edge calculated by Eq. 3.

  * •

When the Kerr medium is aligned at the center of the concave mirrors, KLM can
be only achieved when the intracavity peak power is much lower than the
critical power calculated by Eq. 2.

  * •

To obtain the intracavity power higher than the critical power, the Kerr
medium should be displaced longer than the effective thickness ($d/n_{0}$)
from the center of the concave mirrors.

##  4 Conclusion

We demonstrated the ABCD-matrix analysis of beam propagation in the ring
resonator having the peak power close to or much larger than the critical
power, $P_{\mathrm{c}}$. Our analysis considered the nonlinear propagation
caused by self-focusing in the high-peak-power regime. First, we evaluated the
beam spot size as a function of the resonator parameters, $z$ and $x$ in Fig.
1, and the intracavity power, $P$. When the intracavity power was in the order
of $P_{\mathrm{c}}$, we found a significant dependence on the medium position,
$x$. Moreover, we found that the stable region was continuously deformed with
increasing the intracavity power. The deformation of the stable region
indicated that initially stable CW propagation could be evolved into unstable
pulse propagation and vice versa, depending on the alignment of the concave
mirrors and the Kerr-medium position. Second, we evaluated the pulse stability
by using the $\delta_{K}$ parameter expressed by Eq. 7 to find the optimal
condition of the resonator from the low- to the high-intracavity-power regime.
These two analyses showed the essential fact that the pulses with the peak
power larger than $P_{\mathrm{c}}$ do not favor the condition of the Kerr
medium close to the center of the concave mirrors. On the other hand, the
propagation of these high-peak-power pulses becomes stable by displacing the
Kerr medium from the center position. Therefore, we conclude that stable KLM
can be achieved even for ring resonators with $P/P_{\mathrm{c}}>1$, provided
that they are properly aligned according to the above-mentioned guidelines.

## Acknowledgement

This work was supported by JSPS KAKENHI grant numbers 25707024, 26247049, and
JP17H06205, by the Photon Frontier Network Program of the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), MEXT Quantum Leap
Flagship Program (MEXT Q-LEAP) Grant Number JPMXS0118067246. P. X. was
supported by Advanced Leading Graduate Course for Photon Science (ALPS). We
thank Kazunori Naganuma for his help with the numerical analysis.

## Reference

## References

  * [1] D. E. Spence, P. N. Kean, and W. Sibbett, Opt. Lett. 16, 42–44 (1991).
  * [2] J. M. Evans, D. E. Spence, and W. Sibbett and B. H. T. Chai and A. Miller, Opt. Lett. 17, 1447–1449 (1992).
  * [3] A. Sennaroglu, C. R. Pollock, and H. Nathel, Opt. Lett. 19, 390–392 (1994).
  * [4] V. Petrov, U. Griebner, D. Ehrt, and W. Seeber, Opt. Lett. 22, 408–410 (1997).
  * [5] M. N. Cizmeciyan, H. Cankaya, A. Kurt, and A. Sennaroglu, Opt. Lett. 34, 3056–3058 (2009).
  * [6] T. Brabec, Ch. Spielmann, P. F. Curley, and F. Krausz, Opt. Lett. 17, 1292–1294 (1992).
  * [7] V. Magni, G. Cerullo, and S. D. Silvestri, Opt. Commun. 96, 348–355 (1993).
  * [8] V. Magni, G. Cerullo, and S. D. Silvestri, Opt. Commun. 101, 365–370 (1993).
  * [9] G. Cerullo, S. D. Silvestri, V. Magni, and L. Pallaro, Opt. Lett. 19, 807–809 (1994).
  * [10] V. Magni, J. Opt. Soc. Am. B 13, 2498–2507 (1996).
  * [11] D. Georgiev, J. Herrmann, and U. Stamm, Opt. Commun. 92, 368–375 (1992).
  * [12] T. Brabec, P. F. Curley, Ch. Spielmann, E. Wintner, and A. J. Schmidt, J. Opt. Soc. Am. B 10, 1029–1034 (1993).
  * [13] D. R. Heatley, A. M. Dunlop, and W. J. Firth, Opt. Lett. 18, 170–172 (1993).
  * [14] A. Agnesi, IEEE J. Quantum Electron. 4, 1115–1121 (1994).
  * [15] C. Jirauschek, F. X. Kärtner, and U. Morgner, J. Opt. Soc. Am. B 20, 1356–1368 (2003).
  * [16] L. G. Wright, P. Sidorenko, H. Pourbeyram, Z. M. Ziegler, A. Isichenko, B. A. Malomed, C. R. Menyuk, D. N. Christodoulides, and F. W. Wise, Nature Phys. (2020).
  * [17] V. I. Donin, D. V. Yakovin, and A. V. Gribanov, Opt. Lett. 37, 338–340 (2012).
  * [18] S. Yefet and A. Pe’er, Appl. Sci. 3, 694–724 (2013).
  * [19] A. Klenner and U. Keller, Opt. Express 23, 8532–8544 (2015).
  * [20] Z. Li, J. Peng, J. Yao, and M. Han, Opt. Laser Technol. 23, 1–5 (2017).
  * [21] N. Sugiyama, H. Tanaka, and F. Kannari, Jpn. J. Appl. Phys. 57, 052701 (2018).
  * [22] G. Fibich and A. L. Gaeta, Opt. Lett. 25, 335–337 (2000).
  * [23] M. A. Larotonda, Opt. Commun. B 228, 381–388 (2003).
  * [24] A. Bartels and H. Kurz, Opt. Lett. 27, 1839–1841 (2002).
  * [25] A. Bartels, T. Dekorsy, and H. Kurz, Opt. Lett. 24, 996–998 (1999).
  * [26] T. Brabec, Ch. Spielmann, and F. Krausz, Opt. Lett. 24, 1961–1963 (1991).
  * [27] F. X. Kärtner, J. A. der Au, and U. Keller, IEEE J. Sel. Top. Quantum Electron. 24, 159–168 (1998).
  * [28] H. A. Haus, IEEE J. Sel. Top. Quantum Electron. 6, 1173–1185 (2000).