[1]Dominik Walter<EMAIL_ADDRESS>
# Terahertz microresonators for material characterisation

Dominik Walter Vogt  The Dodd-Walls Centre for Photonic and Quantum
Technologies, Department of Physics, The University of Auckland, Auckland
1010, New Zealand

###### Abstract

Terahertz (THz) technology is rapidly evolving, and the advancement of data
and information processing devices is essential. Silicon THz microresonators
provide perfect platforms to develop compact and integrated devices that could
transform THz technology. Here we present a systematic study on the key figure
of merit of silicon THz disc microresonators - the quality factor (Q-factor) -
in dependence on the substrate’s resistivity. Our results show that the
Q-factor depends linearly on the resistivity and a variation in resistivity
from 10 k$\Omega$cm to 15 k$\Omega$cm changes the Q-factor from 50k to 76k at
0.6 THz. Moreover, we experimentally determine that the silicon material
absorption is inversely proportional to the substrate’s resistivity. In
general, the presented methodology is ideally suited to precisely measure the
material absorption of low-loss materials in the THz domain, which is
challenging using conventional THz spectroscopy.

###### keywords:

Terahertz spectroscopy; microresonator;

## Introduction

Silicon THz microresonators continue to attract significant interest as
fundamental building blocks for essential devices like sensors, filters,
isolators, and temporal differentiators Vogt et al., 2020b ; Yuan et al.,
2019a ; Yuan et al., 2019b ; Gandhi et al., (2021); Zhou et al., (2021); Wang
et al., 2019a ; Yuan et al., (2021). In particular, silicon THz
microresonators stand out due to very high Q-factors of up to 120,000 as well
as the convenient implementation pathway for integrated devices Vogt et al.,
2020a ; Wang et al., 2019b . Various designs ranging from sub-wavelength thin
discs achieving record Q-factors and sensitivities to racetrack ring
resonators ideal for on-chip integration have been explored Vogt et al., 2020a
; Wang et al., 2019b . For example, implementing voltage actuated thermal
tuning and the THz magneto-optical effect in InSb allowed for the recent
realisation of integrated broadly tunable filters and isolators, respectively
Wang et al., 2019a ; Yuan et al., (2021). While significantly more research is
required, these rapid advances highlight the vast potential of silicon THz
microresonator devices in the THz frequency range.

Pivotal to the success of silicon THz microresonators is the low material
absorption of silicon, as the intrinsic losses of the substrate ultimately
limit the Q-factor. The extraordinary low material absorption of high
resistivity float zone grown silicon (HRFZ-Si) in the THz domain has been
known for a long time Xie et al., (2021); Dai et al., (2004); Li and Li,
(2008); in fact, HRFZ-Si is considered lossless for most applications.
However, because of the diminutive absorption it is challenging to precisely
measure the material absorption with conventional methods like THz time-domain
spectroscopy. This is highlighted by the extensive range of results reported
in the literature Xie et al., (2021); Dai et al., (2004); Li and Li, (2008)
and references within. Moreover, while small variations of the silicon
material absorption might be insignificant for many applications, for silicon
THz microresonator devices the material absorption profoundly affects the
Q-factor and thus the device performance.

In this work, we systematically investigate the Q-factor of HRFZ-Si THz disc
microresonators for a range of resistivities and determine the material
absorption of silicon as a function of its resistivity. Our results highlight
how critical resistivity is to achieve high-quality silicon THz microresonator
devices, and allow a straightforward estimate of the achievable Q-factor based
on the substrate’s resistivity. In general, the presented methodology can be
applied to precisely determine the material absorption of low loss materials
in the THz domain.

## Methods

We characterise the Q-factor and resistivity of six HRFZ-Si disc
microresonators fabricated from three wafers (two disc microresonators per
wafer) at 0.5 THz and 0.6 THz, and the experimental results are supported with
finite-element simulations as well as an analytical model. The experiments
performed at 0.5 THz and 0.6 THz serve as a proof-of-concept demonstration;
the presented methodology is applicable to the entire THz domain.

Starting with 2 inch to 4 inch diameter HRFZ-Si wafers, the disc
microresonators are fabricated using femto-second laser micromachining with
220 fs pulses at 800 nm with a repetition rate of 1 kHz and 20 mW power.
Subsequently, the disc’s rim is polished with a fine diamond slurry to remove
imperfections caused by the laser micromachining. The diamond slurry polishing
leads to slightly rounded edges of the discs; however, this and any radial
symmetric features do not deteriorate the intrinsic Q-factors of the
resonators. Moreover, to validate the reproducibility of the fabrication
process, we characterise two discs from each wafer. Figure 1 (a) and (b) show
a typical disc rim before and after the polishing step, respectively. The
polishing is necessary to avoid surface scattering losses which would reduce
the Q-factor, albeit not completely suppress the resonances. Please note that
only the disc’s rim is polished as the laser micromachining does not introduce
imperfections to the top and bottom surfaces of the disc.

Figure 1: (a) 362 $\upmu$m thick disc after the laser micromachining without
polishing. (b) The same disc after several iterative polishing steps with
diamond slurry. Please note these images are not taken at the same position on
the disc’s rim. Both images are taking with an optical microscope with a 20x
magnification objective; the scale bar shows 100 $\upmu$m.

The Q-factor of the THz disc microresonators are characterised using coherent
THz frequency-domain spectroscopy, and data analysis based on Hilbert
transform Vogt and Leonhardt, (2018); Vogt et al., (2019). A schematic of the
experimental setup is shown in Fig. 2 (a). At the heart of the setup is the
single-mode sub-wavelength air-silica step-index waveguide (200 $\upmu$m
diameter) used to evanescently couple the linearly polarised THz radiation to
the disc microresonators. The position of the microresonator relative to the
waveguide is controlled with a 2D computerised translation stage. The typical
coupling distance between the waveguide and disc microresonator depending on
the Q-factor and resonance frequency is about 100 $\upmu$m to 150 $\upmu$m. In
order to avoid absorption from water vapour, the microresonator is placed
inside a gas cell with less than 20 part-per-million water vapour as closely
monitored with a commercial high-end hygrometer. A photograph of a
microresonator close to the waveguide is shown in Fig. 2 (b). The optical
surface quality of the disc’s rim is clearly visible.

Figure 2: (a) Schematic of the experimental setup to study the Q-factor of the
investigated THz disc microresonators. The frequency-domain THz spectrometer
with photo-conductive antennas (PCAs) is a TeraScan 1550nm from Toptica
Photonics Deninger et al., (2015). The high numerical aperture symmetric pass
polymer lenses are specifically designed to coupled to the 200 $\upmu$m
diameter air-silica step index waveguide Vogt and Leonhardt, (2018). (b)
Photograph of a 395 $\upmu$m thick silicon disc in the setup next to the
coupling waveguide; the discs are glued to a thin rod ($\approx$500 $\upmu$m
diameter). The black arrow indicates the direction of propagation of the THz
field while the white arrow indicates the polarisation. (c) A simulated
typical electric field distribution of a higher order radial mode in a 6 mm
diameter 395 $\upmu$m thick THz disc microresonator at 0.6 THz, like the one
shown in (b). Because the mode is confined to the outer 1.5 mm of the disc
microresonator, only 2 mm from the rim of the disc is shown, as indicated with
the black box in (b). The polarisation is parallel to the disc’s rotational
axis.

The intrinsic Q-factor $Q_{i}$ (Q-factor of the uncoupled resonator) is
extracted by fitting an analytical model to the investigated resonances. This
process is repeated for six measurements with slightly different coupling
positions i.e. loaded Q-factors $Q_{l}$, and the average of the extracted
intrinsic Q-factor from those measurements is given in the manuscript. Whereby
the loaded Q-factor $Q_{l}$ is given by
$Q_{l}={({Q_{i}}^{-1}+{Q_{c}}^{-1})}^{-1}$, with $Q_{c}$ the coupling Q-factor
Gorodetsky and Ilchenko, (1999). The detailed description of the THz frequency
domain spectrometer and the data analysis based on the Hilbert transform is
described in great detail elsewhere Vogt et al., (2019).

The resistivity of the discs is measured using a commercial four-point probe
setup with tungsten carbide tips to penetrate the native silica layer of
silicon. The tips are carefully placed around the centre of the disc, and the
finite diameter and thickness of the disc are accounted for in the
calculations Smits, (1958). Resistivity measurements with positive and
negative polarity are performed to ensure ohmic contact with the sample, and
the given resistivities are the averages of the two readings. Interestingly,
the resistivity measurements reveal large variations even on one wafer as
discussed in detail below. For example, a HRFZ-Si wafer specified with a
resisitivity ¿10 k$\Omega$cm has resistivities of 12 k$\Omega$cm and 18
k$\Omega$cm on different positions on the wafer.

The Q-factor measurements are supported with finite element simulations using
COMSOL Multiphysics® software COMSOL AB, (2022). In particular, we model the
fabricated discs using the measured diameter and thickness and adjust the
material absorption until the simulated Q-factor matches the measurements. By
doing so, we can link the measured Q-factor and resistivity to the material
absorption of silicon; this provides a supplementary method to the analytical
model discussed below. The rounded edges from the polishing are not considered
in the model; however, as mentioned above, the rounded edges do not impact the
Q-factor. Please note that the simulations are Eigenfrequency simulations of
the rotational symmetric disc microrseonators in two dimensions and do not
contain the coupling waveguide. Figure 2 (c) shows a typical mode profile of a
higher-order radial mode experimentally excited in the HRFZ-Si THz disc
microresonators. As nearly all of the energy of the different modes is
confined within the discs, it is irrelevant which modes exactly are
experimentally excited, as variations of the Q-factor for different modes are
insignificant compared to the experimentally observed uncertainty.

The finite element simulations also guide the design of the disc
microresonators. In particular, two considerations determine the diameter of
the discs. On the one hand, the diameter has to be large enough so that the
radiation losses are significantly smaller than the losses due to the material
absorption. On the other hand, the diameter should be as small as possible for
a large free-spectral range to avoid difficulties with an over-population of
modes in a narrow frequency range. The intrinsic Q-factor is given as
$Q_{i}={({Q_{m}}^{-1}+{Q_{r}}^{-1})}^{-1}$, with $Q_{m}$ the Q-factor
determined from the material absorption and $Q_{r}$ the radiation Quality
factor. As long as $Q_{r}$ is significantly larger than $Q_{m}$, the intrinsic
Q-factor is determined by the material absorption which is essential for the
presented method. Based on the simulations, a good balance is achieved with 6
mm diameter HRFZ-Si THz disc microresonators with a minimum disc thickness of
200 $\upmu$m. For example, for a 6 mm diameter, and 395 $\upmu$m thick disc
radiation losses are neglectable with a $Q_{r}$ of $6.4\times{10}^{8}$, which
is significantly larger then the maximum anticipated $Q_{m}$ of about
$1\times{10}^{5}$. The disc thicknesses used in this work are 362 $\upmu$m,
395 $\upmu$m, and 523 $\upmu$m - meeting the requirements of neglectable
radiation losses.

Figure 3: (a) Transmission of the air-silica step-index waveguide coupled to a
395 $\upmu$m thick disc microresonator in the frequency range from 592 GHz to
608 GHz, normalised to the waveguide transmission without microresonator. (b)
The corresponding phase profile to (a). (c) High resolution scan of the
resonance at 599 GHz [highlighted in red in (a) and (b)]. (d) Shows the
corresponding phase profile to (c). Please note that the normalised
transmission and phase profile shown in (c) and (d), respectively, are
recorded independently of (a) and (b) at a different coupling position. (e)
and (f) show the normalised transmission and phrase profile of a resonance of
the same THz disc microresonator at 509 GHz. The lower Q-factor compared to
the resonance at 599 GHz is nicely visible, due to the comparable coupling
strengths. The frequency step size in sub-figures (c) to (f) is 1 MHz.

## Results and Discussion

A typical normalised transmission and phase profile of the air-silica step-
index waveguide coupled to one of the 395 $\upmu$m thick disc microresonators
is shown in Fig. 3 (a) and (b) in the frequency range from 592 GHz to 608 GHz,
respectively; a multitude of modes are excited with various coupling strengths
due to different phase matching and evanescent field overlap Vogt and
Leonhardt, (2018). The highlighted mode (in red) at 599 GHz is used to extract
the intrinsic Q-factor of this disc microresonator, and a high resolution
zoom-in is shown in Fig. 3 (c) and (d). As described above, the intrinsic
Q-factor is extracted by fitting an analytical model (solid orange line),
which agrees very well with the measurements Gorodetsky and Ilchenko, (1999).

The measured Q-factors at 0.5 THz (blue crosses) and 0.6 THz (green crosses)
for all discs are summarised in Fig. 4 (a) as a function of the material
absorption. The material absorption is calculated from the measured Q-factor
based on an analytical model Gorodetsky and Ilchenko, (1999):

$Q=\frac{2\pi{n}_{s}f}{\mathrm{c}\alpha},$ (1)

with ${n}_{s}$ the substrate refractive index of 3.416 Vogt et al., (2018), c
the speed of light and $\alpha$ the material absorption. The Q-factor is
inversely proportional to the material absorption, and the corresponding
curves are plotted in solid blue and green lines for 0.5 THz and 0.6 THz,
respectively. We confirmed the validity of the analytical model with finite-
element simulations, where we adjust the material absorption until the
simulated Q-factor matches the experimental observation for the particular
disc microresonator. The analytical model and the simulations [blue and green
circles in Fig. 4 (a)] are almost in perfect agreement. However, the
analytical model is more convenient to use and provides a better
understanding.

Next, Fig. 4 (b) depicts the material absorption as a function of the measured
resistivity $\rho$. As previously suggested, the silicon material absorption
is inversely proportional to the resistivity because the absorption is
predominantly caused by free-charge carriers in silicon Dai et al., (2004).
The experimental results are fitted with $\alpha=b/\rho$, with $b$ the only
fitting parameter, with $b=84.5\,\Omega$ at 0.5 THz and $b=67.5\,\Omega$ at
0.6 THz. The corresponding fits are shown with solid blue and green lines in
Fig. 4 (b). Generally, both frequencies at 0.5 THz (blue crosses) and 0.6 THz
(green crosses) follow the same trend, albeit slightly higher material
absorption at 0.5 THz. This is in agreement with previously reported results
Vogt and Leonhardt, (2018). Moreover, the material absorption shown in Fig. 4
(b) is in agreement with a previously reported upper limit of 0.1 ${cm}^{-1}$
in the frequency range from 0.2 THz to 1 THz Dai et al., (2004). A direct
comparison of the material absorption is unfortunately not possible as
commonly used techniques like THz time-domain spectroscopy can only measure
upper limits for the minute material absorption of HRFZ-Si.

Finally, Fig. 4 (c) shows the measured Q-factors as a function of the
resistivity measured with the four-point probe setup. Because of the inverse
proportionality of the material absorption with resistivity, the Q-factor of
the THz disc microresonators is linear with resistivity (see Eq. 1). The solid
blue and green lines in Fig. 4 (c) are plotted based on the analytical model
and the fitted material absorption $\alpha$ as a function of the resistivity
$\rho$ with fitting parameter $b$.

Figure 4: (a) Q-factor of the THz HRFZ-Si disc microresonators as a function
of the silicon material absorption at 0.5 THz (blue crosses) and 0.6 THz
(green crosses). The corresponding analytical model and the finite-element
simulations are shown with solid lines and circles, respectively. Please note
that two data points for each frequency are very close together and therefore
only five data points per frequency are visible at the shown scale. (b)
Silicon material absorption as a function of the measured resistivity [0.5 THz
(blue) and 0.6 THz (green)]. The functional dependence of the fitted solid
lines is: $\alpha=b/\rho$, with $b=84.5\,\Omega$ at 0.5 THz and
$b=67.5\,\Omega$ at 0.6 THz (c) Q-factor as a function of resistivity at 0.5
THz (blue) and 0.6 THz (green). The solid lines are plotted with Eq. 1 but
using the fit from Fig. 4 (b) for the equation.

Importantly, disc microresonators with similar resistivities have almost
matching Q-factors, confirming the reproducibility of the fabrication process.
Interestingly, however, the 362 $\upmu$m thick discs, while cut from the same
wafer, have resistivities of 12 k$\Omega$cm and 18 k$\Omega$cm, which is
reflected in the measured Q-factors of 86k and 112k, respectively. The
inconsistent resistivities across one wafer are most likely related to minute
amount of imperfections in the wafer itself at or after the fabrication. This
shows, that individual resistivity measurements of commercially available
wafers with a specified resistivity of $>$10 k$\Omega$cm are necessary to
ensure consistent and high-end device performance.

## Conclusions

The presented results highlight the importance of the silicon substrate’s
resistivity on the intrinsic Q-factor of THz microresonators - the building
blocks for next-generation THz devices. Consistent high-end performance
requires rigorous testing of the silicon substrate’s resistivity, as the
slightest variations - even within one wafer - can lead to significant
Q-factor degradation. In particular, our results show that Q-factors of more
than 100k can be achieved by rigorous resistivity measurements of the silicon
substrate.

Moreover, the presented methodology is also suitable to measure the material
absorption of low-loss materials in the THz domain, which can be challenging
using conventional THz spectroscopy. To the best of our knowledge, this is the
first time this technique has been applied in the THz domain to precisely
measure the material absorption of HRFZ-Si as a function of the substrate’s
resistivity.

## Conflicts of interest

The authors declare no conflict of interest.

## References

  * COMSOL AB, (2022) COMSOL AB, Stockholm, S. (2022).  Comsol multiphysics® software.  www.comsol.com.
  * Dai et al., (2004) Dai, J., Zhang, J., Zhang, W., and Grischkowsky, D. (2004).  Terahertz time-domain spectroscopy characterization of the far-infrared absorption and index of refraction of high-resistivity, float-zone silicon.  JOSA B, 21(7):1379–1386.
  * Deninger et al., (2015) Deninger, A. J., Roggenbuck, A., Schindler, S., and Preu, S. (2015).  2.75 THz tuning with a triple-DFB laser system at 1550 nm and InGaAs photomixers.  J. Infrared Millim. Terahertz Waves, 36(3):269–277.
  * Gandhi et al., (2021) Gandhi, R., Leonhardt, R., and Vogt, D. W. (2021).  Microresonator frequency reference for terahertz precision sensing and metrology.  IEEE Trans. Terahertz Sci. Technol., 12(1):70–74.
  * Gorodetsky and Ilchenko, (1999) Gorodetsky, M. L. and Ilchenko, V. S. (1999).  Optical microsphere resonators: optimal coupling to high-Q whispering-gallery modes.  JOSA B, 16(1):147–154.
  * Li and Li, (2008) Li, J. and Li, J. (2008).  Dielectric properties of silicon in terahertz wave region.  Microwave and optical technology letters, 50(5):1143–1146.
  * Smits, (1958) Smits, F. (1958).  Measurement of sheet resistivities with the four-point probe.  Bell System Technical Journal, 37(3):711–718.
  * Vogt et al., (2019) Vogt, D. W., Erkintalo, M., and Leonhardt, R. (2019).  Coherent continuous wave terahertz spectroscopy using Hilbert transform.  Journal of Infrared, Millimeter, and Terahertz Waves, 40(5):524–534.
  * (9) Vogt, D. W., Jones, A. H., Haase, T. A., and Leonhardt, R. (2020a).  Subwavelength thick ultrahigh-Q terahertz disc microresonators.  Photonics Research, 8(7):1183–1188.
  * Vogt et al., (2018) Vogt, D. W., Jones, A. H., and Leonhardt, R. (2018).  Thermal tuning of silicon terahertz whispering-gallery mode resonators.  Applied Physics Letters, 113(1):011101.
  * (11) Vogt, D. W., Jones, A. H., and Leonhardt, R. (2020b).  Terahertz gas-phase spectroscopy using a sub-wavelength thick ultrahigh-Q microresonator.  Sensors, 20(10):3005.
  * Vogt and Leonhardt, (2018) Vogt, D. W. and Leonhardt, R. (2018).  Ultra-high Q terahertz whispering-gallery modes in a silicon resonator.  APL Photonics, 3(5):051702.
  * (13) Wang, Z., Dong, G., Yuan, S., Chen, L., Wu, X., and Zhang, X. (2019a).  Voltage-actuated thermally tunable on-chip terahertz filters based on a whispering gallery mode resonator.  Optics Letters, 44(19):4670–4673.
  * (14) Wang, Z., Yuan, S., Dong, G., Wang, R., Chen, L., Wu, X., and Zhang, X. (2019b).  On-chip single-mode high-Q terahertz whispering gallery mode resonator.  Optics Letters, 44(11):2835–2838.
  * Xie et al., (2021) Xie, J., Ye, W., Zhou, L., Guo, X., Zang, X., Chen, L., and Zhu, Y. (2021).  A review on terahertz technologies accelerated by silicon photonics.  Nanomaterials, 11(7):1646.
  * Yuan et al., (2021) Yuan, S., Chen, L., Wang, Z., Deng, W., Hou, Z., Zhang, C., Yu, Y., Wu, X., and Zhang, X. (2021).  On-chip terahertz isolator with ultrahigh isolation ratios.  Nature communications, 12(1):1–8.
  * (17) Yuan, S., Chen, L., Wang, Z., Wang, R., Wu, X., and Zhang, X. (2019a).  Mode coupling in a terahertz multi-mode whispering-gallery-mode resonator.  Optics letters, 44(8):2020–2023.
  * (18) Yuan, S., Chen, L., Wang, Z., Wang, R., Wu, X., and Zhang, X. (2019b).  Tunable high-quality fano resonance in coupled terahertz whispering-gallery-mode resonators.  Applied Physics Letters, 115(20):201102.
  * Zhou et al., (2021) Zhou, S., Tu, X., Wu, Q., Liu, M., Zhang, Y., Chen, P., Wang, W., Xu, Y., Deng, J., Zhao, Q., et al. (2021).  Photonics-inspired terahertz whispering gallery mode resonator waveguide on silicon platform.  Applied Physics Letters, 119(17):171103.


# Towards Conditional Generation of Minimal Action Potential Pathways for
Molecular Dynamics

John Kevin Cava
&John Vant
&Nicholas Ho
&Ankita Shukla
&Pavan Turaga
&Ross Maciejewski
&Abhishek Singharoy
&
Arizona State University
{jcava,jvant, nichola2, ashukl20, pturaga, rmacieje<EMAIL_ADDRESS>
###### Abstract

In this paper, we utilized generative models, and reformulate it for problems
in molecular dynamics (MD) simulation, by introducing an MD potential energy
component to our generative model. By incorporating potential energy as
calculated from TorchMD into a conditional generative framework, we attempt to
construct a low-potential energy route of transformation between the helix
$\rightarrow$ coil structures of a protein. We show how to add an additional
loss function to conditional generative models, motivated by potential energy
of molecular configurations, and also present an optimization technique for
such an augmented loss function. Our results show the benefit of this
additional loss term on synthesizing realistic molecular trajectories.

##  1 Introduction

The structure and dynamics of proteins are central to our knowledge of
molecular biology. Protein conformational changes are intimately connected to
biological function. For example, some motor proteins changes their local
shape from a helix to a coil-like form to force catalytic sites to change
conformations in binding reacting and releasing sequences. Complementing our
knowledge of 3-dimensional atomic structures, molecular dynamics or MD
simulations has allowed researchers to investigate protein function by
modeling such conformational changes. The two most important parameters for
tracking a conformational change are: a) the transition rate between shapes,
and b) each shape’s relative equilibrium population. Determining these kinetic
and thermodynamic parameters with MD enables a mechanistic description of
protein functions.

Large-scale conformational rearrangements involve non-linear and concerted
movement of many atoms. Recent research in MD simulations has resulted in many
advanced sampling techniques to extract kinetic and thermodynamic parameters
from molecular data. However, in order to extract useful dynamics information,
knowledge of the minimum-action or low-energy pathways is essential. The
energetic barriers crossed by a non-optimal path will cause an erroneous
estimation of the transition rate between states. Some notable examples of
transition path finding algorithms include anisotropic network models [1], the
steepest decent dynamics [2], steered MD [9, 15] and the string simulations
[14]. The two main disadvantages to current path optimization algorithms are:
(i) the computational expense, and, (ii) the complexity and skill required to
set up the calculations. Here, we present a pathway optimization algorithm to
derive the minimum energy pathway between two protein conformational states by
starting from a set of trial pathways.

We utilize a conditional GAN architecture, in which the time-step of molecular
dynamics is used as a conditional for generation. The generator creates a
protein model of vector-size $3N$ for an $N$-atomic system, each atom
described with a feature-size of 3. The potential energy along the backbone
angles of these models are computed using the TorchMD package to constrain the
generation, and also as an input for the discriminator. In order to increase
quality of the model, we performed a pre-training of the generator with
minimization on the dihedral angles and potential energy.

The stretching of deca-alanine from a coiled alpha-helix to an uncoiled
stretched conformational state is used as a model system for protein
conformational changes, due to its small size and relatively simple dynamics.
The generated pathways were compared to ones derived from string with steered
MD simulations and were in good agreement. The github link to the code and
data: https://github.com/johncava/Molecular_Dyanmics_cGAN

##  2 Related Work

A number of approaches are available for learning and generation physical
trajectories using deep learning models. Many previous works that use deep
learning to model physical systems rely heavily on using a model to
approximate a differential equation. Neural ordinary differential equations
use the adjoint sensitivity method to update neural network parameters to
approximate a continuous differential equation [3]. Hamiltonian Neural
Networks (HNNs) predict a scalar value and take the symplectic gradient in
order to approximate a differential equation that conserves total energy [7].
Symplectic-RNN attempts to improve on the stability of HNNs by introducing
recurrent training [4]. Symplectic-ODE Net is a different approach that uses
multiple neural networks to approximate a scalar value to which they take the
symplectic gradient and use adjoint sensitivity to update the parameters [20].
All of these models are useful for modeling non-noisy deterministic few-body
systems, however, MD simulations are inherently stochastic many-body systems
due to thermal fluctuations that follow the Langevin equations. The slow
changes in reaction coordinate (such as collective angle or distances) along a
pathway takes place over a much longer period of time than that of the
atomistic thermal fluctuations. Another approach is to use a deep model to
generate a new frame given a set of history frames using, by using temporal
dynamical models such as LSTMs [18, 8]. Other works have used Langevin
dynamics to find stable conformations [17]. Lastly, methods have used HNNs
with GNNs in order to train on data to enforce a potential energy bias in
order to propel a string of beads in an MD simulation towards a desired shape
(by demonstrating differentiable control) [19]. However, the determination of
the potential energy from the data may not provide a pathway, or even if it is
guaranteed to be minimal. In addition, we want to generalize the model to
predict outside the stable conformation, such that we can hopefully generate
non-equilibrium trajectories.

##  3 Preliminaries

###  3.1 Analytical Potentials

In this paper, we utilize TorchMD [6] to calculate the potential energy terms
of a structure at a given MD time-step. We then use these potentials in
conjunction with cGANs for the generator to learn realistic, low-energy
structural models.

#### Bonds.

The potential energy of the bonds depends on the overall distances between all
the atoms that have an immediate connection to one another. $k_{0}$ is defined
as the force constant and $r$ is the distance between bonded atoms. $r_{eq}$
is the equilibrium distance between the bonded atoms. Then,
$V_{bonded}=k_{0}(r-r_{eq})^{2}$.

#### Angles.

Angles are the second-most important potential energy of the moelcule, in
which torchmd calculates $\theta$ as the angle between the three bonded atoms.
$k_{0}$ is the angular force constant, and $\theta_{eq}$ is the equilibrium
angle. Then, $V_{angle}=k_{\theta}(\theta-\theta_{eq})^{2}$.

#### Dihedrals.

The third important is the dihedral angles, but in TorchMD, it is called
torsion. $\psi$ is the dihedral angle between the four atoms. $\gamma$ is the
phase offset, and $k_{n}$ is the amplitude of the harmonic component of
periodicity $n$. Then,
$V_{torsion}=\sum_{n=1}^{n_{max}}{k_{n}(1+cos(n\psi-\gamma))}$.

#### Lennard Jones.

In this paper, we consider this Lennard Jones, whereas TorchMD calls this Van
der Waals energy. In previous work, such as Ham-Net [12], they use this
equation as a loss for their model, but where A and B is proportional to the
output weights of their deep learning model. Then,
$V_{VdW}=\frac{A}{r^{12}}-\frac{B}{r^{6}}$, where, $A=4e\sigma^{12}$ and
$B=4e\sigma^{2}$.

#### Electrostatics.

Electrostatics is the Coulomb’s force, such that $k_{e}=\frac{1}{4\pi e_{0}}$,
and $q_{i}$ and $q_{r}$ are the charges of two atoms. Then,
$V_{electrostatics}=k_{0}\frac{q_{i}q_{r}}{r}$, where, $A=4e\sigma^{1}2$ and
$B=4e\sigma^{2}$.

#### Total Energy.

TorchMD analytically calculates all the previous potential energies
individually, and outputs a dictionary in which we can add all the potential
energy to become the total potential energy for a given structure. We also
added an external potential energy (End2End distance) that adds a force bias
(given a timestep t) that we want the generated trajectory to follow. Overall,
the total potential energy is defined as
$V_{total}=\sum^{n_{bonds}}{V_{bonded}}+\sum^{n_{angles}}{V_{angles}}+\sum^{n_{torsion}}{V_{torsion}}\\\
+\sum^{n_{atoms}}_{i}\sum^{n_{atoms}}_{j<i}{(V_{VdW}+V_{electrostatics})}+V_{End2End}(t)$.

##  4 Experiments and Results

Our experimental implementation is available in Pytorch [16] and TorchMD,
which was executed on a computing cluster with NVIDIA Tesla K80s.

#### Data and Data Augmentation

Figure 1: cGAN Model. A Generator that takes in a time step and random vector,
which outputs a frame that is then used as an input to the Discriminator and
the calculation of the potential energy.

For our experiments, we utilize the data for simulating a 10 deca-alanine
molecule given a bias force that would stretch it. For these experiments, we
conducted an augmentation such that it is randomly rotated and all the
trajectories are aligned on the center of mass of the molecule. We did this in
order to make sure that the generator does not over-fit on specific positions
and translations.

#### Conditional GAN with Geometric Invariances.

We take inspiration from MolGAN [5], cGANs [13], and InvNet [10] in order to
generate structures that are conditioned on time. The loss function is
augmented by the potential energy computed through TorchMD. Our model is
described in Figure 1. The generator is a feed forward neural network that
takes in a time step, and a normal vector of size 31. This then generates the
structure which would be 40 atoms (in backbone prediction) and 104 atoms (in
full system). We also add a constraint by computing the potential energy from
the generated structures in order to train a better generator. The generator
is a feedforward neural network: [Linear(32,50), Linear(50,100),
Linear(100,312)]. The discriminator takes 1 time-step + 250 torsion-angles
computed from the 312 sized output structure. Thus, the discriminator is
structured as another feedforward neural network: [Linear(251,10),
Linear(10,1)]. For the optimizer, we used Adam [11] for optimizing the
generator and discriminator (learning rate = $1\times 10^{-3}$). The number of
epochs for GAN training is 10.

Figure 2: Four generated frames trained by cGAN alongside the RMSD plots with
respect to the ground truth trajectory.

#### Loss Function.

The loss function we use is a standard cGAN loss, augmented with the total
potential energy term described in section 3.1. The loss function is then,

$\min_{G}\max_{D}L(D,G)+\lambda V_{total},$ (1)

where, $L(D,G)$ is the standard cGAN loss function, conditioned on time-step
$t$, given by $\mathbb{E}_{x\sim p_{data}(x)}[logD(x|t)]+\mathbb{E}_{x\sim
p_{z}(z)}[log[1-D(G(z|t))]$. Due to the additional term $V_{total}$, which
actually depends on the generator parameters, a total derivative of the loss
function is difficult. Motivated by the three-way update rule that was used in
InvNet [10], we find that a similar strategy works well even in our case. That
is, a GAN-like update of $G$ via gradient steps of $L(D,G)$ keeping $D$ fixed;
a GAN-like update of $D$ via gradient steps of $L(D,G)$ keeping $G$ fixed; and
an update of $G$ via gradient steps of $V_{total}$.

#### Pre-Training the GANs with Annealing.

In order to narrow down the search space, we implement pretraining of the
generator before we do the adversarial training. In addition, we conduct an
annealing process for the pretraining. In this process, we slowly transition
from training on reconstruction (of position/angles) to potential energies
(bonds, angles, dihedrals). In the reconstruction loss, we have an emphasis on
torsion angles such that we emphasize the initial structure of the helix for
the deca-alanine. Thus once we do the adversarial training, we then emphasize
on the End2End potential energy bias force exerted on the deca-alanine. The
number of epochs for pretraining is 5, and we used a learning rate of $1\times
10^{-2}$.

#### Adversarial Training and Potential Encoding.

In our adversarial training, we emphasize an invariant conditional component
(as inspired by physical constraint GANs such as invGAN). Instead of a neural
network, we use potential energies calculated from TorchMD given the generated
structures. We have experimented with a hierarchy of potential energies given
a certain epoch. Where we focus in order of bonds, angles, dihedrals,
impropers, LJ, and electrostatics. Together, they enforce the physical
constraints that make the generator learn how to make a physically realistic
structure. We have also made it such that with pretraining we have the
generator creating good enough structure based on the bonds, angles, and
dihedrals, since LJ is sensitive to plainly wrong structures, which would
produce high potential energy, and thus blow up the gradients. In addition, to
avoid blowing up of gradients, we did gradient clipping on the generator such
that during training the potential energy loss wouldn’t blow up the gradients
for the generator.

#### Results and Conclusion.

After training the cGAN, we generated 20 time points conditioned on the
maximum trajectory size of 1002 frames. In Figure 2, we present some exemplary
generated time points from our cGAN model. The lowest point of each RMSD plot
is the ‘ground truth match’ to the generated frame. We find that as the time
points progress from 4 to 9 to 14, our generated models match frames 6, 8 and
11 of the ground truth data (downsampled from 1002 to 20 frames). Based on
these frames, we can conclude that the model is able to find the minimial
action potential pathway from a favorable conformational state to a less
conformational state. We have tried training the cGAN without the $V_{total}$,
and we were not able to generate physically realistic structures.

However, there are still limitations in which the starting structure at time
point 0 is not helical, which indicates that there needs to be better ways of
representing structure as adding a discriminator and potential energy loss
constraint aren’t enough to enforce that initial structure. Future work could
be considering representing the input as a ball and stick model and treat
changes as ‘actions’ in a reinforcement learning framework, or also
representing the input as a spatio-temporal graph.

## References

  * [1] Ali Rana Atilgan, SR Durell, Robert L Jernigan, Melik C Demirel, O Keskin, and Ivet Bahar.  Anisotropy of fluctuation dynamics of proteins with an elastic network model.  Biophysical journal, 80(1):505–515, 2001.
  * [2] Changjun Chen, Yanzhao Huang, Xiaofeng Ji, and Yi Xiao.  Efficiently finding the minimum free energy path from steepest descent path.  The Journal of chemical physics, 138(16):164122, 2013.
  * [3] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  arXiv preprint arXiv:1806.07366, 2018.
  * [4] Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou.  Symplectic recurrent neural networks.  arXiv preprint arXiv:1909.13334, 2019.
  * [5] Nicola De Cao and Thomas Kipf.  Molgan: An implicit generative model for small molecular graphs.  arXiv preprint arXiv:1805.11973, 2018.
  * [6] Stefan Doerr, Maciej Majewsk, Adrià Pérez, Andreas Krämer, Cecilia Clementi, Frank Noe, Toni Giorgino, and Gianni De Fabritiis.  Torchmd: A deep learning framework for molecular simulations, 2020.
  * [7] Sam Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  CoRR, abs/1906.01563, 2019.
  * [8] Chitrak Gupta, John Kevin Cava, Daipayan Sarkar, Eric A Wilson, John Vant, Steven Murray, Abhishek Singharoy, and Shubhra Kanti Karmaker.  Mind reading of the proteins: Deep-learning to forecast molecular dynamics.  bioRxiv, 2020.
  * [9] Sergei Izrailev, Sergey Stepaniants, Barry Isralewitz, Dorina Kosztin, Hui Lu, Ferenc Molnar, Willy Wriggers, and Klaus Schulten.  Steered molecular dynamics.  In Computational molecular dynamics: challenges, methods, ideas, pages 39–65. Springer, 1999.
  * [10] Ameya Joshi, Minsu Cho, Viraj Shah, Balaji Pokuri, Soumik Sarkar, Baskar Ganapathysubramanian, and Chinmay Hegde.  Invnet: Encoding geometric and statistical invariances in deep generative models.  Proceedings of the AAAI Conference on Artificial Intelligence, 34(04):4377–4384, Apr. 2020.
  * [11] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [12] Ziyao Li, Shuwen Yang, Guojie Song, and Lingsheng Cai.  Conformation-guided molecular representation with hamiltonian neural networks.  In International Conference on Learning Representations, 2020.
  * [13] Mehdi Mirza and Simon Osindero.  Conditional generative adversarial nets.  CoRR, abs/1411.1784, 2014.
  * [14] Albert C Pan, Deniz Sezer, and Benoît Roux.  Finding transition pathways using the string method with swarms of trajectories.  The journal of physical chemistry B, 112(11):3432–3440, 2008.
  * [15] Sanghyun Park and Klaus Schulten.  Calculating potentials of mean force from steered molecular dynamics simulations.  The Journal of chemical physics, 120(13):5946–5961, 2004.
  * [16] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  CoRR, abs/1912.01703, 2019.
  * [17] Chence Shi, Shitong Luo, Minkai Xu, and Jian Tang.  Learning gradient fields for molecular conformation generation.  arXiv preprint arXiv:2105.03902, 2021.
  * [18] Sun-Ting Tsai, En-Jui Kuo, and Pratyush Tiwary.  Learning molecular dynamics with simple language model built upon long short-term memory neural network.  Nature communications, 11(1):1–11, 2020.
  * [19] Wujie Wang, Simon Axelrod, and Rafael Gómez-Bombarelli.  Differentiable molecular simulations for control and learning, 2020.
  * [20] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty.  Symplectic ode-net: Learning hamiltonian dynamics with control.  arXiv preprint arXiv:1909.12077, 2019.


$\|\bm{\Delta}_{\delta}\|_{*}\chi(\|\bm{\Delta}_{\delta}\|_{*}\leq
N^{-\frac{3}{K}})\lesssim\left(|\Theta_{\delta}|+N^{\frac{2}{K}}\|\mathbf{D}\|_{*}\right)\chi(\|\bm{\Delta}_{\delta}\|_{*}\leq
N^{-\frac{3}{K}})\prec N^{-\frac{5}{K}}.$ (5.42)

We have thus found a “forbidden” area for $\|\bm{\Delta}_{\delta}\|_{*}$. With
the aid of a standard stochastic continuity argument, we remove the indicator
function from the bounds (5.41) and (5.42) to get to the rough bounds

$\|\bm{\Delta}_{\delta}\|_{*}\prec
N^{-\frac{5}{K}}\quad\text{and}\quad|\Theta_{\delta}|\prec N^{-\frac{5}{K}}.$
(5.43)

Since $\mathbf{x}$ and $\mathbf{y}$ were arbitrary we the first bound becomes
$|\bm{\Delta}_{\delta}|\prec N^{-5/K}.$ Now, assume that
$\max\\{|\bm{\Delta}_{0}|,|\bm{\Delta}_{1}|\\}\prec\Lambda$ and
$\max\\{|\Theta_{0}|,|\Theta_{1}|\\}\prec\theta$ for some deterministic
$\theta\leq\Lambda\leq N^{-3/K}$. Then we know from Lemma 5.5 and Proposition
4.5 as well as (5.17) that

$\max_{i}\\{|\bm{\Delta}_{i}|\\}\prec\theta+N^{\frac{2}{K}}\sqrt{\frac{\operatorname{\mathrm{Im}}m+\Lambda}{N\eta}}\quad\text{and}\quad\max_{i}\\{|\xi_{1}\Theta_{i}+\xi_{2}\Theta_{i}^{2}|\\}\prec
N^{\frac{2}{K}}\frac{\operatorname{\mathrm{Im}}m+\Lambda}{N\eta}$ (5.44)

are also deterministic bounds. Here, we also used
$\operatorname{\mathrm{Im}}m\sim\langle\operatorname{\mathrm{Im}}M_{1}\rangle$
by Corollary 3.3, Lemma 4.1 as well as (4.18). The first bound in (5.44) self-
improving and applying it iteratively yields

$\max_{i}\\{|\bm{\Delta}_{i}|\\}\prec\theta+N^{\frac{2}{K}}\left(\frac{N^{\frac{2}{K}}}{N\eta}+\sqrt{\frac{\operatorname{\mathrm{Im}}m+\theta}{N\eta}}\right)$
(5.45)

and therefore the second bound in (5.44) improves to

$\max_{i}\\{|\xi_{1}\Theta_{i}+\xi_{2}\Theta_{i}^{2}|\\}\prec
N^{\frac{2}{K}}\frac{\operatorname{\mathrm{Im}}m+\theta}{N\eta}+N^{\frac{4}{K}}\frac{1}{(N\eta)^{2}}.$
(5.46)

Now we separately treat
$\operatorname{\mathrm{Re}}z\in\operatorname{supp}(\rho)$ and
$\operatorname{\mathrm{Re}}z\notin\operatorname{supp}(\rho)$ and we start with
the former. Then by Corollary 3.3 we know that
$\operatorname{\mathrm{Im}}m\sim\sqrt{\kappa+\eta}$. For fixed $\theta$ we
apply Lemma 5.6 with

$d=N^{\frac{2}{K}}\frac{\sqrt{\kappa+\eta}+\theta}{N\eta}+N^{\frac{4}{K}}\frac{1}{(N\eta)^{2}}$
(5.47)

to obtain

$\max_{i}\\{|\Theta_{i}|\\}\prec\min\left\\{\frac{d}{\sqrt{\kappa+\eta}},\sqrt{d}\right\\}.$
(5.48)

This is also a self-improving bound and iterating it gives

$\max_{i}\\{|\Theta_{i}|\\}\prec
N^{\frac{2}{K}}\frac{1}{N\eta},\quad\text{hence}\quad\max_{i}\\{\bm{|}\Delta_{i}|\\}\prec
N^{\frac{2}{K}}\left\\{\sqrt{\frac{\operatorname{\mathrm{Im}}m}{N\eta}}+\frac{N^{\frac{2}{K}}}{N\eta}\right\\}.$
(5.49)

For $\operatorname{\mathrm{Re}}z\notin\operatorname{supp}(\rho)$ we have
$\operatorname{\mathrm{Im}}m\sim\frac{\eta}{\sqrt{\kappa+\eta}}$, again by
Corollary 3.3. Analogously to the
$\operatorname{\mathrm{Re}}z\in\operatorname{supp}(\rho)$ we obtain

$\max_{i}\\{|\Theta_{i}|\\}\prec
N^{\frac{2}{K}}\frac{1}{N(\eta+\kappa)}+N^{\frac{4}{K}}\frac{1}{(N\eta)^{2}\sqrt{\kappa+\eta}}.$
(5.50)

Finally we use (5.20), (4.31) and the bounds on $\max_{i}\\{|\Theta_{i}|\\}$
to arrive at the averaged bounds

$\max_{i}\\{|\bm{\Delta}_{i}|_{\mathrm{av}}\\}\prec
N^{\frac{2}{K}}\begin{cases}\frac{1}{N\eta}\text{ if
}\operatorname{\mathrm{Re}}z\in\mathrm{supp}{\rho_{1}}\\\
\frac{1}{N(\eta+\kappa)}+\frac{N^{\frac{2}{K}}}{N^{2}\eta^{2}(\eta+\kappa)^{\frac{1}{2}}}\text{
if }\operatorname{\mathrm{Re}}z\notin\mathrm{supp}{\rho_{1}}.\end{cases}$
(5.51)

This concludes the proof of (5.27) on $\mathbb{D}_{\gamma_{1}}^{\kappa_{0}}$.

Now, assume (5.28) on $\mathbb{G}_{\gamma_{0}}^{C,\eta_{0}}$. From (4.33) we
have

$\mathscr{L}[\mathbf{\Delta}_{\delta}]=-M_{\delta}\mathbf{D}+M_{\delta}\mathcal{S}[\mathbf{\Delta}_{\delta}]\mathbf{\Delta}_{\delta}.$
(5.52)

We apply $\mathscr{L}^{-1}$ to the equation and take its *-norm for some
deterministic $\mathbf{x},$ $\mathbf{y}$ to find

$\begin{split}\|\mathbf{\Delta}_{\delta}\|_{*}&\leq\|\mathscr{L}^{-1}[M_{\delta}\mathbf{D}]\|_{*}+\|\mathscr{L}^{-1}[M_{\delta}\mathcal{S}[\mathbf{\Delta}_{\delta}]\mathbf{\Delta}_{\delta}]\|_{*}\\\
&\leq\|\mathscr{L}^{-1}\|_{*\to*}\left(\|M_{\delta}\mathbf{D}\|_{*}+\|M_{\delta}\mathcal{S}[\mathbf{\Delta}_{\delta}]\mathbf{\Delta}_{\delta}\|_{*}\right)\\\
&\lesssim_{C}N^{\frac{2}{K}}\|\mathbf{D}\|_{*}+N^{\frac{2}{K}}\|\mathbf{\Delta}_{\delta}\|_{*}^{2}.\end{split}$
(5.53)

Here, $\|\mathscr{L}^{-1}\|_{*\to*}$ denotes the operator norm of
$\mathscr{L}^{-1}$ with respect to the *-norm. In the last estimate, we have
used
$\|\mathscr{L}^{-1}\|_{*\to*}\lesssim\|\mathscr{L}^{-1}\|_{\mathrm{sp}}\lesssim_{C}1$
by [22, Equation (70c)] and Proposition 4.6 as well as
$\|M_{\delta}\mathbf{R}\|_{*}\lesssim N^{2/K}\|\mathbf{R}\|_{*}$ and
$\|M_{\delta}\mathcal{S}[\mathbf{R}]\mathbf{R}\|_{*}\lesssim
N^{2/K}\|\mathbf{R}\|_{*}^{2}$ by [22, Equation (70a) and (70b)]. Equation
(5.42) also holds on $\mathbb{G}_{\gamma_{1}}^{C,\eta_{0}}$ by the same
argument as in the case $z\in\mathbb{D}_{\gamma_{1}}^{\kappa_{0}}$ and for
sufficiently large $K$ and sufficiently small $\gamma_{s}$ we have in
particular $\|\mathbf{D}\|_{*}\prec N^{-\frac{7}{K}}$ for
$\delta\in\\{0,1\\}$. We use this bound on with (5.53) to estimate

$\|\mathbf{\Delta}_{\delta}\|_{*}\chi(\|\mathbf{\Delta}_{\delta}\|_{*}\leq
N^{-\frac{3}{K}})\prec_{C}N^{-\frac{5}{K}}.$ (5.54)

By a stochastic continuity argument, we establish the rough bound

$\|\mathbf{\Delta}_{\delta}\|_{*}\prec_{C}N^{-\frac{5}{K}}\quad\text{and
therefore}\quad|\mathbf{\Delta}_{\delta}|\prec_{C}N^{-\frac{5}{K}}$ (5.55)

since $\mathbf{x}$ and $\mathbf{y}$ were arbitrary. Now, assume
$\max_{i}\\{|\bm{\Delta}_{i}|\\}\prec_{C}\Lambda$ for some deterministic
$\Lambda\leq N^{-3/K}$. Then we have another deterministic bound in the form
of

$\max_{i}\\{|\bm{\Delta}_{i}|\\}\prec_{C}N^{\frac{2}{K}}\sqrt{\frac{\operatorname{\mathrm{Im}}m+\Lambda}{N\eta}}.$
(5.56)

Iterating this self-improving bound and using
$\operatorname{\mathrm{Im}}m\sim_{C}\eta$ from Lemma A.3 we find

$\max_{i}\\{|\bm{\Delta}_{i}|\\}\prec_{C}\frac{N^{\frac{4}{K}}}{N\eta}+N^{\frac{2}{K}}\frac{1}{\sqrt{N}}.$
(5.57)

Finally from (5.57), (5.52) and (4.31) we obtain the averaged bound

$\max_{i}\\{|\mathbf{\Delta}_{i}|_{\mathrm{av}}\\}\prec_{C}\frac{N^{\frac{4}{K}}}{(N\eta)^{2}}+\frac{1}{N}.$
(5.58)

Therefore we have shown that (5.28) holds on
$\mathbb{G}_{\gamma_{1}}^{C,\eta_{0}}$. ∎

###### Proof of Propositions 4.3 and 4.4.

For all $K>0$ Equation (5.27) holds on $\mathbb{D}_{1}^{\kappa_{0}}$ for some
$\kappa_{0}$ and (5.28) on $\mathbb{G}_{1}^{C,\eta_{0}}$ for all $C>0$ and
some $\eta_{0}$ depending on $C$ due to the global law Proposition 5.7. The
local laws, Proposition 4.3 and Proposition 4.4, follow immediately from
applying Proposition 5.8 finitely many times in the respective domains. ∎

This concludes the proof of Theorem 2.9 and Proposition 2.10 and we are left
with proving their Corollaries, 2.20 and 2.12.

###### Proof of Corollary 2.20.

Let $\lambda_{i}$, $i\in\llbracket N\rrbracket$, denote the eigenvalues of
$q(\mathbf{X})$ with eigenvectors $\mathbf{v}_{i}$. let
$\mathbf{x}\in\mathbb{C}^{N\times N}$ be a deterministic and normalized vector
and let $\mathbf{v}$ be an eigenvector of $q(\mathbf{X})$ with eigenvalue
$\lambda$ such that $|\lambda-\tau_{0}|<\kappa_{0}$ for $\kappa_{0}$ from
Theorem 2.9. Evaluating $\mathbf{g}$ at spectral parameter
$\lambda+\mathrm{i}\eta$ have with very high probability that

$1\gtrsim\operatorname{\mathrm{Im}}\langle\mathbf{x},\mathbf{g}(\lambda+\mathrm{i}\eta)\mathbf{x}\rangle=\sum_{i=1}^{N}\eta\frac{\eta}{\eta^{2}+(\lambda-\lambda_{i})^{2}}|\langle\mathbf{x},\mathbf{v}_{i}\rangle|^{2}\geq\frac{|\langle\mathbf{x},\mathbf{v}_{i}\rangle|}{\eta}$
(5.59)

for all $\eta\geq N^{-1+\varepsilon}$ and all $\varepsilon>0$. As the
deterministic vector $\mathbf{x}$ was arbitrary, the eigenvalue
delocalization, (2.20), follows. ∎

###### Proof of Corollary 2.12.

Let $q$ have regular edge at $\tau_{0}$. W.l.o.g. we can assume that
$\tau_{0}=\tau_{+}$ is a right edge. Otherwise, consider the right edge of
$-q$. From (2.18) and Proposition 2.10 it follows that

$\mathbb{P}\left(|\langle\mathbf{B}(\mathbf{g}-m)\rangle|>\|\mathbf{B}\|N^{\varepsilon}\left(\frac{1}{N(\kappa+\eta)}+\frac{1}{(N\eta)^{2}\sqrt{\kappa+\eta}}\right)\right)\lesssim_{\varepsilon,\gamma,D,C}N^{-D}$
(5.60)

for all deterministic $\mathbf{B}\in\mathbb{C}^{N\times N}$ and
$z\in\mathbb{D}_{\gamma}^{\kappa_{0}}\cup\mathbb{G}_{\gamma}^{C,\eta_{0}}$
with $E\notin\operatorname{supp}(\rho)$. In particular, this holds true for
$C\geq\kappa_{0}^{-1}$. Following [24, Chapter 11.1], this implies for all
$\varepsilon>0$ that with very high probability there are no eigenvalues
$\lambda$ such that $\lambda\notin\operatorname{supp}(\rho)$ and
$N^{-2/3+\varepsilon}\leq\lambda-\tau_{+}\leq C$, i.e.

$\mathbb{P}\left(\exists\lambda\in\operatorname{Spec}(q(\mathbf{X}))\,:\,N^{-\frac{2}{3}+\varepsilon}\leq\lambda-\tau_{+}<C\right)\lesssim_{\varepsilon,D,C}N^{-D}.$
(5.61)

Additionally, it follows from the trivial bound (4.32) that there is a $K>0$
such that

$\mathbb{P}\left(\|q(\mathbf{X})\|\geq K\right)\lesssim_{D}N^{-D}.$ (5.62)

Combining (5.61) for some $C>0$ such that $C+\tau_{0}\geq K$ and (5.62) we
have

$\mathbb{P}\left(\exists\lambda\in\operatorname{Spec}(q(\mathbf{X}))\,:\,\lambda-\tau_{+}\geq
N^{-\frac{2}{3}+\varepsilon}\right)\lesssim_{\varepsilon,D}N^{-D}.$ (5.63)

By a standard argument, (5.63) in conjunction with the averaged local law,
(2.17), implies eigenvalue rigidity around the edge as in Theorem 2.12 (see
e.g. [24, Chapter 11.2-11.4]). ∎

##  Appendix A Appendix

###  A.1 The entrywise real part of a Hermitian matrix.

Let $n\in\mathbb{N}$ and $H\in\mathbb{C}^{n\times n}$ be a Hermitian matrix
and $\widehat{H}=\frac{1}{2}(H+H^{t})$ be its entrywise real part. The
following lemma summarizes how the two matrices are related.

######  Lemma A.1.

Let $h_{i}$ and $\widehat{h}_{i}$ be the eigenvalues of $H$ and $\widehat{H}$
respectively, both arranged in non-increasing order. The following holds true

  1. 1.

$h_{1}\geq\widehat{h}_{1}$ and $h_{n}\leq\widehat{h}_{n}$.

  2. 2.

$\|H\|\geq\|\widehat{H}\|$.

  3. 3.

If $H\geq 0$, then $\widehat{H}\geq 0$ and conversely $H\leq 0$ implies
$\widehat{H}\leq 0$.

  4. 4.

Let additionally $\operatorname{rank}H=1$. Then we have
$\operatorname{rank}\widehat{H}=1$ if and only if $H\in\mathbb{R}^{l\times
l}$. Otherwise we have $\operatorname{rank}\widehat{H}=2$.

###### Proof.

Let $\widetilde{H}=\frac{1}{2\mathrm{i}}(H-H^{t})$ be the entrywise imaginary
part of $H,$ i.e. $H=\widehat{H}+\mathrm{i}\widetilde{H}$. First, we will
prove $h_{1}\geq\widehat{h}_{1}$. Note that for all $a\in\mathbb{R}$ the norms
of $H+a$ and $\widehat{H}+a$ are given by $\|H+a\|=\max\\{h_{1}+a,-h_{n}-a\\}$
and $\|\widehat{H}+a\|=\max\\{\widehat{h}_{1}+a,-\widehat{h}_{n}-a\\}$. Thus
for sufficiently large $a\in\mathbb{R}$ (depending on $H$) we have
$\|H+a\|=h_{1}+a$ and $\|\widehat{H}+a\|=\widehat{h}_{1}+a$. For such an $a$
let $v_{1}\in\mathbb{R}^{n}$ be a normalized eigenvector of $\widehat{H}+a$
corresponding to the eigenvalue $\widehat{h}_{1}+a$ (it can be chosen purely
real since $\widehat{H}$ is a real symmetric matrix). Then

$h_{1}+a=\|H+a\|\geq\|(H+a)v_{1}\|=\sqrt{\|(\widehat{H}+a)v_{1}\|^{2}+\|\widetilde{H}v_{1}\|^{2}}\geq\|(\widehat{H}+a)v_{1}\|=\widehat{h}_{1}+a.$
(A.1)

Here the second equality holds since $\widetilde{H}v_{1}$ is a purely
imaginary vector. The claim $h_{1}\geq\widehat{h}_{1}$ follows. Choosing $a$
sufficiently small, we find $h_{n}\leq\widehat{h}_{n}$ by a similar argument.
Since $\|H\|=\max\\{h_{1},-h_{n}\\}$ and
$\|\widehat{H}\|=\max\\{\widehat{h}_{1},-\widehat{h}_{n}\\}$, the inequality
$\|H\|\geq\|\widehat{H}\|$ follows. Next, let $H\geq 0$. Then $h_{n}\geq 0$
and thus $\widehat{h}_{n}\geq h_{n}\geq 0$. The inequality $\widehat{H}\geq 0$
follows. Similarly, $H\leq 0$ implies $\widehat{H}\leq 0$.

Now let $\operatorname{rank}H=1$. Then $\operatorname{rank}\widehat{H}=1$ if
$H\in\mathbb{R}^{n\times n}$ is clear. Let on the other hand
$H\in\mathbb{C}^{n\times n}\setminus\mathbb{R}^{n\times n}$. Then $H=\alpha
vv^{*}$ for some $\alpha\in\mathbb{R}\setminus\\{0\\}$ and normalized
$v\in\mathbb{C}^{n}\setminus e^{\mathrm{i}\varphi}\mathbb{R}^{n}$ for all
$\varphi\in\mathbb{R}$. Then
$\widehat{H}=\frac{\alpha}{2}(vv^{*}+\bar{v}\bar{v}^{*})$ and since $v$ and
$\bar{v}$ are only linearly dependent if $v\in
e^{\mathrm{i}\varphi}\mathbb{R}^{n}$ the claim
$\operatorname{rank}\widehat{H}=2$ follows. ∎

###  A.2 Matrix inequality

Here we provide a proof for the matrix inequality used in (5.15).

######  Lemma A.2.

Let $R,T\in\mathbb{C}^{n\times n}$ be arbitrary matrices. Then the following
inequality holds:

$(R+T)^{*}(R+T)\leq 2(R^{*}R+T^{*}T).$ (A.2)

###### Proof.

We first note that

$(R+T)^{*}(R+T)=R^{*}R+T^{*}T+R^{*}T+T^{*}R$ (A.3)

and it is thus sufficient to bound the last two terms. Let
$v\in\mathbb{C}^{n}$ be arbitrary. Since $R^{*}T+T^{*}R$ is Hermitian,
$\langle v,(R^{*}T+T^{*}R)v\rangle\in\mathbb{R}$ and we estimate

$\langle v,(R^{*}T+T^{*}R)v\rangle=2\operatorname{\mathrm{Re}}(\langle
Rv,Tv\rangle)\leq 2\|Rv\|\|Tv\|\leq\|Rv\|^{2}+\|Tv\|^{2}=\langle
v,(R^{*}R+T^{*}T)v\rangle,$ (A.4)

where we used the Cauchy-Schwarz inequality in the second step. As $v$ was
arbitrary, $R^{*}T+R^{*}T\leq R^{*}R+T^{*}T$ follows and we conclude the
lemma. ∎

###  A.3 Some properties of the solution of the Dyson Equation

###### Proof of Proposition 2.3.

There can at most be one analytic function in the upper half-plane that
satisfies (2.4) and $\lim_{z\to\infty}zm(z)=-1$ since (2.4) is stable at
infinity. Thus we are left with proving existence. First, consider the case of
$A$ being invertible. Let $\\{s_{1},\ldots,s_{l}\\}$ be a family of free semi-
circular variables in a $C^{*}$ probability space $(\mathscr{S},\tau)$ (see
[21, Appendix B]) and let $s:=(s_{i})_{i\in\llbracket l\rrbracket}$. Define

$\mathbf{L}_{\mathrm{sc}}=K_{0}\otimes\mathbb{1}+\sum_{j=1}^{l}K_{j}\otimes
s_{j}.$ (A.5)

Following the proof of [21, Lemma 2.6], the matrix-valued function
$M_{0}:\mathbb{H}\to\mathbb{C}^{(l+1)\times(l+1)},$ given by

$M_{0}(z)=(\mathrm{id}\otimes\tau)(\mathbf{L}_{\mathrm{sc}}-zJ\otimes\mathbb{1})^{-1},$
(A.6)

is a solution to (4.14) for $\delta=0$. Thus its (1,1) entry satisfies (4.20)
for $\delta=0$ and therefore (2.4). Using the Schur complement formula, we
find

$(M_{0}(z))_{11}=\tau(q(s)-z\otimes\mathbb{1})^{-1}.$ (A.7)

The polynomial $q(s)$, defined in (1.1) is self-adjoint, thus $(M_{0})_{11}$
is analytic on $\mathbb{H}$, has positive imaginary part and
$\lim_{z\to\infty}z(M_{0}(z))_{11}=-1$. Therefore $m:=(M_{0})_{11}$ is the
unique function that satisfies all conditions of Proposition 2.3.

Now, let $A$ be non-invertible. Then there is a $u\sim 1$ such that
$A^{\varepsilon}:=A+\varepsilon$ is invertible for all $0<\varepsilon<u.$ Let
$q^{\varepsilon}$ and $\gamma^{\varepsilon}$ be the objects given by replacing
$A$ with $A+\varepsilon$ in the definitions of $q$ in (1.1) and $\gamma$ in
(2.5), respectively. By the above argument, the function

$m^{\varepsilon}(z):=\tau(q^{\varepsilon}(s)-z\otimes\mathbb{1})^{-1}$ (A.8)

is a solution to the equation

$-\frac{1}{m^{\varepsilon}}=z+\gamma^{\varepsilon}(m^{\varepsilon}).$ (A.9)

For any fixed $z\in\mathbb{H}$, both $m^{\varepsilon}$ and
$\gamma^{\varepsilon}$ are continuous in $\varepsilon$ at $\varepsilon=0$.
Thus $m:=m^{0}$ solves (2.4) and $m$ is analytic on $\mathbb{H}$, has positive
imaginary part and satisfies $\lim_{z\to\infty}zm(z)=-1$. Therefore it is the
unique function that satisfies all conditions of Proposition 2.3. ∎

The function $m$ is a Stieltjes transform of a real-valued probability measure
$\rho$ with compact support. As such it can be analytically extended to
$\mathbb{C}\setminus\operatorname{supp}(\rho).$ The following lemma summarizes
the properties of the extension, also called $m$, on and above the real axis
outside of the spectrum

######  Lemma A.3.

For the analytical extension of $m$ to
$\mathbb{C}\setminus\operatorname{supp}(\rho)$ we have the following.

  1. 1.

The function $m$ is real-valued on
$\mathbb{R}\setminus\operatorname{supp}(\rho)$ and $m^{\prime}(E)>0$ for all
$E\in\mathbb{R}\setminus\operatorname{supp}(\rho)$.

  2. 2.

For all $C>0$ there is an $\eta_{0}>0$ such that for all $z=E+\mathrm{i}\eta$
with $C^{-1}\leq\operatorname{dist}(E,\operatorname{supp}(\rho))\leq C$ and
$0<\eta\leq\eta_{0}$ we have $\operatorname{\mathrm{Im}}m\sim_{C}\eta.$

###### Proof.

The fact that $m(E)\in\mathbb{R}$ for all
$\mathbb{R}\setminus\operatorname{supp}(\rho)$ follows immediately for from
(2.6) and taking the derivative of (2.6) we find

$m^{\prime}(E)=\int_{\mathbb{R}}\frac{\rho(\mathrm{d}x)}{(x-E)^{2}}>0.$ (A.10)

Using $m(E)\in\mathbb{R}$ we have for $z\in\mathbb{H}$

$\operatorname{\mathrm{Im}}m(z)=\operatorname{\mathrm{Im}}(m(z)-m(E))=\eta
m^{\prime}(E)+\mathcal{O}(\eta^{2}).$ (A.11)

By (A.10) and continuity of the derivative, we have that $m^{\prime}(E)$ is
bounded from above and bounded away from 0 on all compact subsets of
$\mathbb{R}\setminus\operatorname{supp}(\rho)$. Therefore we have
$m^{\prime}(E)\sim_{C}1$ for all $C>0$ and $E$ such that
$C^{-1}\leq\operatorname{dist}(E,\operatorname{supp}(\rho))\leq C$. ∎

###### Proof of Lemma 4.1.

The function $m$ is a Stieltjes transform and as such analytically extendable
to $\mathbb{C}\setminus\operatorname{supp}(\rho)$. Since, by assumption,
$\rho$ has a regular edge at $\tau_{\pm}$ it can also be continuously extended
to $\tau_{\pm}$. For
$x\in(\mathbb{R}\setminus\operatorname{supp}(\rho))\cup\\{\tau_{\pm}\\}$ we
denote this extension by $m(x).$ By continuity it satisfies

$-\frac{1}{m(x)}=x+\gamma(m(x)).$ (A.12)

Consider (4.20) at $z\in\mathbb{H}$ as a perturbation to (A.12). More
precisely, consider the following equation in the unknown function
$\widetilde{m}_{\delta}(z)$,

$\frac{1}{m(x)}-\frac{1}{\widetilde{m}_{\delta}(z)}=\gamma_{\delta}(\widetilde{m}_{\delta}(z))-\gamma_{0}(m(x))+(z-x),$
(A.13)

where we have used $\gamma_{0}=\gamma.$ Define by
$\dot{\gamma}_{\delta}(y):=\partial_{t}\gamma_{\frac{t}{\mathrm{i}\eta}}(y)|_{t=\mathrm{i}\delta\eta}$
the partial derivative with respect to $\mathrm{i}\delta\eta$ and by
$\gamma_{\delta}^{\prime}(y):=\partial_{y}\gamma_{\delta}(y)$ the derivative
with respect to its variable.

Expanding (A.13) in $m_{\delta}(z)-m(x)$ and $\mathrm{i}\delta\eta$ we find

$\begin{split}&\left(\frac{1}{m(x)^{2}}-\gamma_{0}^{\prime}(m(x))\right)(\widetilde{m}_{\delta}(z)-m(x))-\left(\frac{1}{m(x)^{3}}+\frac{1}{2}\gamma_{0}^{\prime\prime}(m(x))\right)(\widetilde{m}_{\delta}(z)-m(x))^{2}\\\
&\qquad\qquad=\mathrm{i}\dot{\gamma}_{0}(m(x))\delta\eta+(z-x)+\mathcal{O}(|\widetilde{m}_{\delta}(z)-m(x)|^{3}+|\widetilde{m}_{\delta}(z)-m(x)|\eta+\eta^{2}).\end{split}$
(A.14)

We identify $\frac{1}{m^{2}}-\gamma_{0}^{\prime}(m)=h(m)$ and
$-\frac{2}{m^{3}}-\gamma_{0}^{\prime\prime}(m)=h^{\prime}(m)$ with $h$ defined
in (3.10). Now let $x=\tau_{\pm}$ and note that $m(\tau_{\pm})=m_{\pm}.$ By
Lemma 3.5 as well as Lemma 3.6 we have

$h(m_{\pm})=0,\quad|h^{\prime}(m_{\pm})|\sim 1\quad\text{and}\quad\pm
h^{\prime}(m_{\pm})>0.$ (A.15)

Furthermore one finds $\dot{\gamma}_{0}(y)>0$ for all
$y\in\mathbb{R}\setminus\mathscr{S}(\gamma)$. Thus (A.14) constitutes an
approximate quadratic equation in $\widetilde{m}_{\delta}(z)-m_{\pm}$ and
there are $u,u_{1}\sim 1$ such that for all $z\in\mathbb{H}$ with
$0<|z-\tau_{\pm}|\leq u$, there are exactly two solutions
$\widetilde{m}_{\delta}(z)$ with $|\widetilde{m}_{\delta}(z)-m_{\pm}|\leq
u_{1}$. One of them has positive imaginary part and one of them has negative
imaginary part. Combining (A.14) and (A.15), we find that the unique solution
with positive imaginary part satisfies both parts of (4.22) and we denote it
by $m_{\delta}$.

Now let $C>0$, $x=E=\operatorname{\mathrm{Re}}z$ and
$C^{-1}<\operatorname{dist}(E,\operatorname{supp}(\rho))<C$. Then we have

$h(m(E))=(m^{\prime}(E))^{-1}\sim_{C}1$ (A.16)

by Lemma A.3 and the fact that
$C^{-1}<\operatorname{dist}(E,\operatorname{supp}(\rho))<C$ constitutes a
compact subset of $\mathbb{R}\setminus\operatorname{supp}(\rho)).$ Therefore
(A.14) is an approximate linear equation in $\widetilde{m}_{\delta}(z)-m(E)$.
Hence there are $u,u_{1}\sim 1$ such that for all $0<|z-E|=\eta\leq u$ there
is a unique solution $\widetilde{m}_{\delta}(z)$ with
$|\widetilde{m}_{\delta}(z)-m(E)|\leq u_{1}$ to (A.14). We denote this
solution by $m_{\delta}$ and combining (A.15) with (A.16) we find that it
satisfies both parts of (4.23).

From now on let $A$ be invertible and

$M_{\delta}(z)=(\mathrm{id}\otimes\tau)(\mathbf{L}_{\mathrm{sc}}-(zJ+\delta\eta(I-J))\otimes\mathbb{1})^{-1},$
(A.17)

where $\mathbf{L}_{\mathrm{sc}}$ was defined in (A.5). For all $\eta>0$, the
imaginary part of $zJ+\delta\eta(I-J)$ is positive. Therefore, by Lemma [31,
Lemma 5.4], the function $M_{\delta}$ is the unique solution to (4.14) with
$\operatorname{\mathrm{Im}}M_{\delta}(z)>0$. In particular, this implies that
its (1,1) entry satisfies (4.20). By the Schur complement formula, it is given
by

$(M_{\delta}(z))_{11}=\tau(q_{\delta}(s)-z\otimes\mathbb{1})^{-1},$ (A.18)

where $q_{\delta}$ denotes the polynomial $q$, defined in (1.1), with $A$
replaced by $A_{\delta}$. For $E\notin\operatorname{supp}(\rho)$ consider

$\begin{split}(M_{\delta})_{11}-m&=\tau((q_{\delta}(s)-z\otimes\mathbb{1})^{-1}-(q_{0}(s)-z\otimes\mathbb{1})^{-1})\\\
&=\tau((q_{\delta}(s)-z\otimes\mathbb{1})^{-1}(q_{0}(s)-q_{\delta}(s))(q_{0}(s)-z\otimes\mathbb{1})^{-1}).\end{split}$
(A.19)

The middle factor satisfies

$\|q_{0}(s)-q_{\delta}(s)\|=\left\|\sum_{i,j=1}^{l}(A-A_{\delta})s_{i}s_{j}\right\|\lesssim\eta\delta.$
(A.20)

Since $E\notin\operatorname{supp}(\rho)$, the term $q_{0}(s)-E$ is invertible
and by continuity in $\eta$ we have

$\|(q_{0}(s)-z\otimes\mathbb{1})^{-1}\|\lesssim_{E}1$ (A.21)

uniformly in $\eta$ for $0<\eta\leq 1$. In particular, for $\eta$ sufficiently
small (depending on $E$), we have

$\|(q_{0}(s)-z\otimes\mathbb{1})^{-1}(q_{0}(s)-q_{\delta}(s))\|\leq\frac{1}{2}.$
(A.22)

Therefore the first factor in (A.19) factorizes further into

$(q_{\delta}(s)-z\otimes\mathbb{1})^{-1}=(1\otimes\mathbb{1}+(q_{0}(s)-z\otimes\mathbb{1})^{-1}(q_{0}(s)-q_{\delta}(s)))^{-1}(q_{0}(s)-z\otimes\mathbb{1})^{-1}$
(A.23)

and

$\|(1\otimes\mathbb{1}+(q_{0}(s)-z\otimes\mathbb{1})^{-1}(q_{0}(s)-q_{\delta}(s)))^{-1}\|\leq
2.$ (A.24)

Combining (A.19) with (A.20),(A.21),(A.23) and (A.24), we find

$|(M_{\delta})_{11}-m|\lesssim_{E}\delta\eta.$ (A.25)

As shown in the first part of the proof, there are $u,u_{1}>0$ such that
(A.13) has a unique solution $m_{\delta}$ with positive imaginary part and
$|m_{\delta}(z)-m_{\pm}|\leq u_{1}$ for all $z\in\mathbb{H}$ such that
$|z-\tau_{\pm}|\leq u$. For such $u,u_{1}$ we choose a $u_{2}$ with
$0<u_{2}<u$ such that $|m_{z}-m_{\pm}|\leq\frac{u_{1}}{2}$ for all
$z\in\mathbb{H}$ with $|z-\tau_{\pm}|\leq u_{2}$ (such a $u_{2}$ exists by
Corollary 3.3). Fix $E\notin\operatorname{supp}(\rho)$ with
$|E-\tau_{\pm}|<u_{2}$. We choose $\eta=\operatorname{\mathrm{Im}}z$
sufficiently small such that $|z-\tau_{\pm}|<u_{2}$ and
$|(M_{\delta})_{11}-m|\lesssim\frac{u_{1}}{2}$ (which is possible by (A.25)).
Thus there is a $z\in\mathbb{H}$ with $|z-\tau_{\pm}|\leq u$ such that

$|(M_{\delta})_{11}(z)-m_{\pm}|\leq|(M_{\delta})_{11}(z)-m(z)|+|m(z)-m_{\pm}|\leq
u_{1}.$ (A.26)

Since this condition is unique among solutions to (4.20) with positive
imaginary part, we must have $(M_{\delta})_{11}(z)=m_{\delta}(z)$ for some
$z\in\mathbb{H}$ with $|z-\tau_{\pm}|\leq u$. By continuity of both
$m_{\delta}$ and $(M_{\delta})_{11}$ we must have
$(M_{\delta})_{11}(z)=m_{\delta}(z)$ for all $z\in\mathbb{H}$ with
$|z-\tau_{\pm}|\leq u$.

For $C>0$ and $E$ with
$C^{-1}<\operatorname{dist}(E,\operatorname{supp}(\rho))<C$ we find similarly
from Lemma A.3 and (A.25) that for all $u_{1}>0$ there is an $\eta>0$ such
that

$|(M_{\delta})_{11}(z)-m(E)|\leq|(M_{\delta})_{11}(z)-m(z)|+|m(z)-m(E)|\leq
u_{1}.$ (A.27)

Thus $(M_{\delta})_{11}(z)$ must be the unique solution to (4.20) that
satisfies (4.23), which concludes the proof of the lemma. ∎

## References

  * [1] Johannes Alt, László Erdős, Torben Krüger, and Dominik Schröder.  Correlated random matrices: band rigidity and edge universality.  Ann. Probab., 48(2):963–1001, 2020.
  * [2] Greg W. Anderson.  Convergence of the largest singular value of a polynomial in independent Wigner matrices.  Ann. Probab., 41(3B):2103–2181, 2013.
  * [3] Greg W. Anderson.  A local limit law for the empirical spectral distribution of the anticommutator of independent Wigner matrices.  Annales de l’I.H.P. Probabilités et statistiques, 51(3):809–841, 2015.
  * [4] Greg W. Anderson, Alice Guionnet, and Ofer Zeitouni.  An introduction to random matrices, volume 118 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 2010.
  * [5] Fanny Augeri.  Large deviations principle for the largest eigenvalue of Wigner matrices without Gaussian tails.  Electron. J. Probab., 21:No. 32, 2016.
  * [6] Fanny Augeri, Alice Guionnet, and Jonathan Husson.  Large deviations for the largest eigenvalue of sub-Gaussian matrices.  Comm. Math. Phys., 383(2):997–1050, 2021.
  * [7] Z. D. Bai and Y. Q. Yin.  Necessary and sufficient conditions for almost sure convergence of the largest eigenvalue of a Wigner matrix.  Ann. Probab., 16(4):1729–1741, 1988.
  * [8] Zhidong Bai and Jack W. Silverstein.  Spectral analysis of large dimensional random matrices.  Springer Series in Statistics. Springer, New York, second edition, 2010\.
  * [9] Marwa Banna and Tobias Mai.  Hölder continuity of cumulative distribution functions for noncommutative polynomials under finite free Fisher information.  J. Funct. Anal., 279(8):108710, 2020.
  * [10] Zhigang Bao, László Erdős, and Kevin Schnelli.  Local law of addition of random matrices on optimal scale.  Comm. Math. Phys., 349(3):947–990, 2017.
  * [11] Zhigang Bao, László Erdős, and Kevin Schnelli.  Spectral rigidity for addition of random matrices at the regular edge.  J. Funct. Anal., 279(7):108639, 2020.
  * [12] Serban. T. Belinschi and Mireille Capitaine.  Spectral properties of polynomials in independent Wigner and deterministic matrices.  J. Funct. Anal., 273(12):3901–3963, 2017.
  * [13] Serban T. Belinschi, Tobias Mai, and Roland Speicher.  Analytic subordination theory of operator-valued free additive convolution and the solution of a general random matrix problem.  J. Reine Angew. Math., 732:21–53, 2017.
  * [14] Paul Bourgade, László Erdös, and Horng-Tzer Yau.  Edge universality of beta ensembles.  Comm. Math. Phys., 332(1):261–353, 2014.
  * [15] Ian Charlesworth and Dimitri Shlyakhtenko.  Free entropy dimension and regularity of non-commutative polynomials.  J. Funct. Anal., 271(8):2274–2292, 2016.
  * [16] Benoît Collins, Alice Guionnet, and Félix Parraud.  On the operator norm of non-commutative polynomials in deterministic matrices and iid GUE matrices.  Camb. J. Math., 10(1):195–260, 2022.
  * [17] Nicholas A. Cook, Alice Guionnet, and Jonathan Husson.  Spectrum and pseudospectrum for quadratic polynomials in Ginibre matrices.  Ann. Inst. Henri Poincaré Probab. Stat., 58(4):2284–2320, 2022\.
  * [18] Ken Dykema.  On certain free product factors via an extended matrix model.  J. Funct. Anal., 112(1):31–60, 1993.
  * [19] László Erdős, Antti Knowles, and Horng-Tzer Yau.  Averaging fluctuations in resolvents of random band matrices.  Ann. Henri Poincaré, 14(8):1837–1926, 2013.
  * [20] László Erdős, Antti Knowles, Horng-Tzer Yau, and Jun Yin.  The local semicircle law for a general class of random matrices.  Electron. J. Probab., 18:No. 59, 2013.
  * [21] László Erdős, Torben Krüger, and Yuriy Nemish.  Local laws for polynomials of Wigner matrices.  J. Funct. Anal., 278(12):108507, 2020.
  * [22] László Erdős, Torben Krüger, and Dominik Schröder.  Random matrices with slow correlation decay.  Forum Math. Sigma, 7:No. e8, 2019.
  * [23] László Erdős and Yuanyuan Xu.  Small deviation estimates for the largest eigenvalue of Wigner matrices.  Bernoulli, 29(2):1063–1079, 2023.
  * [24] László Erdős and Horng-Tzer Yau.  A dynamical approach to random matrix theory, volume 28 of Courant Lecture Notes in Mathematics.  Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 2017.
  * [25] László Erdős, Horng-Tzer Yau, and Jun Yin.  Rigidity of eigenvalues of generalized Wigner matrices.  Adv. Math., 229(3):1435–1515, 2012.
  * [26] Ohad N. Feldheim and Sasha Sodin.  A universality result for the smallest eigenvalues of certain sample covariance matrices.  Geom. Funct. Anal., 20(1):88–123, 2010.
  * [27] Fritz Gesztesy and Eduard Tsekanovskii.  On matrix-valued Herglotz functions.  Math. Nachr., 218:61–138, 2000.
  * [28] Friedrich Götze, Alexey Naumov, and Alexander Tikhomirov.  Local laws for non-Hermitian random matrices and their products.  Random Matrices Theory Appl., 9(4):2150004, 2020.
  * [29] Friedrich Götze and Alexander Tikhomirov.  On the Asymptotic Spectrum of Products of Independent Random Matrices.  ArXiv:1012.2710, 2011.
  * [30] Uffe Haagerup, Hanne Schultz, and Steen Thorbjørnsen.  A random matrix approach to the lack of projections in $C^{*}_{\rm red}(\mathbb{F}_{2})$.  Adv. Math., 204(1):1–83, 2006.
  * [31] Uffe Haagerup and Steen Thorbjørnsen.  A new application of random matrices: ${\rm Ext}(C^{*}_{\rm red}(F_{2}))$ is not a group.  Ann. of Math. (2), 162(2):711–775, 2005.
  * [32] J. William Helton, Tobias Mai, and Roland Speicher.  Applications of realizations (aka linearizations) to free probability.  J. Funct. Anal., 274(1):1–79, 2018.
  * [33] J. William Helton, Scott A. McCullough, and Victor Vinnikov.  Noncommutative convexity arises from linear matrix inequalities.  J. Funct. Anal., 240(1):105–191, 2006.
  * [34] J. William Helton, Reza Rashidi Far, and Roland Speicher.  Operator-valued semicircular elements: solving a quadratic matrix equation with positivity constraints.  Int. Math. Res. Not. IMRN, 2007(22):rnm086, 2007.
  * [35] Benjamin Landon and Horng-Tzer Yau.  Edge statistics of dyson brownian motion.  arXiv: 1712.03881, 2017.
  * [36] Ji Oon Lee and Kevin Schnelli.  Edge universality for deformed Wigner matrices.  Rev. Math. Phys., 27(8):1550018, 2015.
  * [37] Ji Oon Lee and Jun Yin.  A necessary and sufficient condition for edge universality of Wigner matrices.  Duke Math. J., 163(1):117–173, 2014.
  * [38] Tobias Mai, Roland Speicher, and Moritz Weber.  Absence of algebraic relations and of zero divisors under the assumption of full non-microstates free entropy dimension.  Adv. Math., 304:1080–1107, 2017.
  * [39] Camille Male.  The norm of polynomials in large random and deterministic matrices.  Probab. Theory Related Fields, 154(3-4):477–532, 2012.  With an appendix by Dimitri Shlyakhtenko.
  * [40] Yuriy Nemish.  Local law for the product of independent non-Hermitian random matrices with independent entries.  Electron. J. Probab., 22:No. 22, 2017.
  * [41] Sean O’Rourke and Alexander Soshnikov.  Products of independent non-hermitian random matrices.  Electron. J. Probab., 16:No. 81, 2219–2245, 2011.
  * [42] Natesh S. Pillai and Jun Yin.  Edge universality of correlation matrices.  Ann. Statist., 40(3):1737–1763, 2012.
  * [43] Natesh S. Pillai and Jun Yin.  Universality of covariance matrices.  Ann. Appl. Probab., 24(3):935–1001, 2014.
  * [44] Dimitri Shlyakhtenko and Paul Skoufranis.  Freely independent random variables with non-atomic distributions.  Trans. Amer. Math. Soc., 367(9):6267–6291, 2015.
  * [45] Alexander Soshnikov.  Universality at the edge of the spectrum in Wigner random matrices.  Comm. Math. Phys., 207(3):697–733, 1999.
  * [46] Terence Tao and Van Vu.  Random matrices: universality of local eigenvalue statistics up to the edge.  Comm. Math. Phys., 298(2):549–572, 2010.
  * [47] Terence Tao and Van Vu.  Random matrices: sharp concentration of eigenvalues.  Random Matrices Theory Appl., 2(3):1350007, 2013.
  * [48] Craig A. Tracy and Harold Widom.  Level-spacing distributions and the Airy kernel.  Phys. Lett. B, 305(1-2):115–118, 1993.
  * [49] Craig A. Tracy and Harold Widom.  On orthogonal and symplectic matrix ensembles.  Comm. Math. Phys., 177(3):727–754, 1996.
  * [50] Dan Voiculescu.  Limit laws for random matrices and free products.  Invent. Math., 104(1):201–220, 1991.
  * [51] Eugene P. Wigner.  Characteristic vectors of bordered matrices with infinite dimensions.  Ann. of Math. (2), 62:548–564, 1955.