Microscopy and Microanalysis

Paper

[$\ast$]Corresponding author<EMAIL_ADDRESS>
# Live Iterative Ptychography

Dieter Weber0000-0001-6635-9567  Simeon Ehrig0000-0002-8218-3116  Andreas
Schropp0000-0000-0000-0000  Alexander Clausen0000-0002-9555-7455  Silvio
Achilles0000-0001-7244-7854  Nico Hoffmann  Michael
Bussmann0000-0002-8258-3881  Rafal E. Dunin-Borkowski0000-0001-8082-0647
Christian G. Schroer0000-0002-9759-1200  Ernst Ruska-Centre for Microscopy and
Spectroscopy with Electrons, Forschungszentrum Jülich, 52425 Jülich, Germany
Center for X-ray and Nano Science CXNS, Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  Helmholtz Imaging, Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany  Helmholtz-Zentrum
Dresden-Rossendorf, D-01328 Dresden, Germany  Center for Advanced Systems
Understanding, Untermarkt 20, 02826 Görlitz, Germany

(2023)

###### Abstract

We demonstrate live-updating ptychographic reconstruction with ePIE, an
iterative ptychography method, during ongoing data acquisition. The
reconstruction starts with a small subset of the total data, and as the
acquisition proceeds the data used for reconstruction is extended. This
creates a live-updating view of object and illumination that allows monitoring
the ongoing experiment and adjusting parameters with quick turn-around. This
is particularly advantageous for long-running acquisitions. We show that such
a gradual reconstruction yields interpretable results already with a small
subset of the data. We show simulated live processing with various scan
patterns, parallelized reconstruction, and real-world live processing at the
hard X-ray ptychographic nanoanalytical microscope PtyNAMi at the PETRA III
beamline.

###### keywords:

ptychography, X-ray microscopy

##  1 Introduction

Ptychography (Hoppe, 1969a; Hoppe & Strube, 1969; Hoppe, 1969b) can
reconstruct a quantitative object transfer function for a specimen from a set
of measurements where a spatially modulated illumination is shifted relative
to the object and the intensity distribution of the transmitted and projected
beam, typically far-field diffraction patterns, are measured for the different
shift positions (Figure 1). The alternating projection method, i.e., PIE/ePIE
reconstruction algorithm (Rodenburg & Faulkner, 2004; Maiden & Rodenburg,
2009) is one of the classical and well-known approaches to perform this
reconstruction.

This method is currently in use at the hard X-ray ptychographic nanoanalytical
microscope PtyNAMi at the PETRA III beamline P06 at DESY (Schropp et al.,
2020). It allows to retrieve images of objects at a higher spatial resolution
than the size of the X-ray focus and is physically only limited by the
collection angle, wave length of the probe, and the available coherent photon
flux (Schropp & Schroer, 2010).

Figure 1: Schematic setup for far field ptychography. The illumination is
localized and scanned over the object O(r). The detector captures the
transmitted beam in the far field for each scan position. Illustration based
on (Bangun et al., 2022).

Here, a ptychographic data set is used for live processing that is recorded in
the standard data format available the beamline. The raw data including the
diffractions patterns and position values are stored in individual hdf5 files
following in general the NeXus format.

Currently, the data is exported from the the standard data format available at
the beamline, and then processed offline using a stand-alone software. The raw
data, including the diffractions patterns and position values, are stored in
individual HDF5 files following, in general, the NeXus format (Könnecke et
al., 2015). (Weber et al., 2023) contains the dataset “scan_00063” with such a
raw beamline data structure. Adjusting the parameters for ptychography
requires a range of test measurements to adjust exposure, beam position on the
detector, specimen region, focal distance and settings for the ptychography
engine. Since many specimens exhibit mostly phase contrast, even simple
previews of the data require a phase contrast method.

Live reconstruction shortens the time until feedback is available to the
operators and can help to use the scarce beam time efficiently by re-adjusting
or aborting measurements before they complete. Measurement times at PtyNAMi
range between minutes for test measurements and many hours for large
acquisitions.

Previously, Strauch et al. showed live single side-band (SSB) ptychography
(Strauch et al., 2021). However, this algorithm is not suitable for this
particular application since the reconstruction resolution is limited to the
scan grid, and it can only reconstruct successfully if the scan step size is
smaller than the resolving power of the illuminating beam, determined by its
convergence angle. The low acquisition speed and resolving power of X-ray
microscopes compared to electron microscopes means reconstruction with higher
resolution than the scan grid, which allows coarse sampling of the specimen,
is highly favorable here. Live Wigner distribution deconvolution (WDD) (Bangun
et al., 2022) has similar limitations.

Near live ptychographic reconstruction has also been shown by Pelz et al.
(Pelz et al., 2021) by speeding up the reconstruction after fast scans so that
it appears nearly instantaneous. However, this is not applicable for this use
case since the acquisition speed and not the reconstruction is the limiting
factor, making a partial reconstruction during an ongoing scan desirable.

Using the well-established extended ptychographical iterative engine (ePIE)
implementation at DESY is favorable since it is in practical use for offline
processing and the users’ know-how is transferable. Towards this goal, three
steps were undertaken:

  1. 1.

Create an API for the established stand-alone ePIE implementation that allows
to decouple data source, ptychographic processing and display of results.

  2. 2.

Investigate the behavior of iterative phase retrieval on partial input data in
a simulated setting to determine if live iterative ptychography is possible
for typical datasets acquired at PtyNAMi.

  3. 3.

Implement live ptychographic reconstruction at PtyNAMi.

##  2 Materials and Methods

A reconstruction with an iterative method encompasses the following aspects:

  * •

Determine general parameters such as object and illumination size in pixel
units, relation between physical coordinates and pixel coordinates.

  * •

Preprocessing in order to match the different coordinate systems of object,
illumination, scan and detector so that the forward and backward propagation
model used in the optimization routine matches the actual measurements
quantitatively.

  * •

Assign initial values to object and illumination.

  * •

Preprocess input data.

  * •

Iterative optimization of illumination and object.

The following subsections detail how these aspects are adapted for live
reconstruction.

###  2.1 ePIE Implementation

For the reconstruction of the experimental data, a GPU accelerated
implementation of the ePIE (Maiden & Rodenburg, 2009) algorithm is used. The
update strength for object and illumination, typically denoted by $\alpha$ and
$\beta$, was 1 for all reconstructions shown here. The software is written in
CUDA C++ and developed by the X-ray Nanoscience and X-ray Optics group at DESY
for PtyNAMi. The ePIE is an iterative algorithm, which reconstructs the object
and the illumination at the same time. The implementation starts with a
general initialization phase which is then followed by the main reconstruction
routine. The initialization phase prepares the diffraction patterns for the
reconstruction and sets up the object and illumination with an initial guess.
In this implementation the diffraction patterns are cropped to a size of
$2^{n}$, with $n$ being an integral value, in order to use a simple fast-
Fourier transformation and increase the performance. Furthermore, the
implementation can be configured to transform the diffraction patterns. Since
the forward and back projection in this ePIE engine is implemented with a
Discrete Fourier Transform (DFT), the sampling and coordinate system of the
reconstructed object transfer function must match the sampling and coordinate
system of the diffraction patterns used in the reconstruction to allow forward
and backward projection with a DFT. In the ePIE implementation used at
PtyNAMi, the spatial resolution and size of the object reconstruction in real
space is determined by the size and scale of the diffraction patterns and the
scan area. That means the spatial resolution and size of the object
reconstruction can be adjusted by cropping the diffraction patterns before
reconstruction. Furthermore, depending on the detector and its orientation,
the recorded raw diffraction patterns might be rotated and/or flipped relative
to the coordinate system of the position encoders. This can be compensated by
rotating or flipping the diffraction patterns in such a way that their
coordinate system and handedness match the position encoders. In the second
phase, the object and illumination are reconstructed. The implementation
iterates over the diffraction patterns in a randomly shuffled order and
calculates an update for the object and illumination. The reconstruction ends
after a specific number of iterations.

###  2.2 Extensions of the ePIE implementation

The original ePIE application is designed to perform an offline
reconstruction, which means all diffraction patterns are available at the
beginning of the reconstruction and the reconstruction configuration remains
constant during runtime. There is no user input to add new diffraction
patterns, and changes of the reconstruction parameters during the
reconstruction are not possible. However, a live reconstruction requires that
diffraction patterns can be added to follow the progressing scan during the
reconstruction. Furthermore, it is necessary to stop and continue the
reconstruction in order to extract intermediate updates, include new data and
match the iteration speed with the acquisition speed to not over-optimize
prematurely on a small subset of the input data.

For the interactive live reconstruction, Python bindings have been added to
the existing C++ reconstruction engine using Pybind11 in order to connect it
with the Python-based data handling and control logic. Compared to C++, Python
provides a read–eval–print loop (REPL) interpreter, which is perfectly suited
for interactive workflows. Furthermore, Jupyter Notebooks can be used, giving
some extra benefits such as an easy programming interface in the web browser
and the possibility to combine source code with different media types, like
markdown text and images to better document the reconstruction.

###  2.3 Preprocessing

Real-world detector data may require preprocessing such as cropping, scaling
and orientation changes. Since the bindings don’t support exporting or
extending the loaded data, but require creating a new instance of the engine
whenever the input data changes, a LiberTEM user-defined function (UDF) was
implemented to perform preprocessing equivalent to the ePIE implementation
based on its configuration file. This allows reusing the preprocessed data and
reduces the amount of data loaded into the ePIE implementation. Furthermore,
LiberTEM can perform this efficiently on live data streams in a parallel
fashion.

###  2.4 Determine parameters

For offline processing of an entire dataset, the ePIE engine determines the
reconstruction parameters from the size of the detector data in pixels after
preprocessing, bounding box of physical scan positions, and optical parameters
such as camera length, detector pixel size and wavelength from the
configuration file. For live processing, the ePIE engine is initialized with
mock detector data and approximate scan positions based on the scan settings
that are known beforehand. The engine uses the nominal scan positions to
calculate the nominal object size and pads it with half the illumination size
so that the full area of the object that intersects with the illumination is
covered. During the experiment, the scan positions are measured by
interferometers or encoders.

Based on this initialization, a configuration for the engine is derived from
the initial configuration that fixes the size and relation between physical
coordinates and pixel coordinates. That way the input data can change
dynamically during live processing without influencing the reconstruction
parameters.

###  2.5 Simulating live and parallel processing

Live processing starts with a small initial subset of data that is extended
gradually as the scan proceeds (Figure 2). As a preliminary test before
implementing true live processing, this was simulated on offline data using
the bindings described above: The initial values for object and illumination
together with a small subset of the input data are loaded, a small number of
iterations are performed, the resulting object and illumination are recorded,
and the cycle is repeated with a larger subset of the input data using the
previous state of object and illumination as initialization values until the
complete dataset is consumed. Different scan schemes such as random sampling
or spiral scans can be simulated by re-ordering the input data and positions.

Figure 2: Schematic progress of live iterative ptychography. Data is
preprocessed and accumulated, and the reconstruction algorithm includes more
and more data in the iteration as it becomes available. This creates an
evolving view of the object and illumination reconstruction. In real-world
live processing, the sequence is determined by the acquisition speed,
iteration speed and possibly limits on the number of iterations to avoid over-
optimization. For simulated live processing, the sequence of available data
and iteration steps is pre-determined by a script.

Parallel processing can be simulated in a similar fashion by independent
iteration on disjoint subsets of the input data with intermittent averaging
between the object and illumination for synchronization, similar to (Nashed et
al., 2014). This creates a hybrid between classical ePIE where the object and
illumination are updated dynamically with each diffraction pattern during a
pass, and batch-oriented algorithms that accumulate updates and apply them
only after a pass.

A ptychographic dataset for offline reconstruction had been recorded at
PtyNAMi (Schropp et al., 2020) on an NTT-AT resolution test chart with 50nm
thinnest lines and spaces (model ATN/XRESO-50HC) at an X-ray energy of E =
13.5 keV. The X-ray beam was focused by a pair of nanofocusing refractive
lenses (NFLs) (Schroer et al., 2005; Schropp et al., 2010) to a spot size of
116 nm (horizontal) x 118 nm (vertical). The sample was positioned slightly in
defocus at a distance of approximately 200um. The diffraction patterns were
collected with an Eiger X 4M detector (Dectris AG) positioned at a distance of
3.2m behind the sample. In the reconstructions these patterns were cropped to
a size of (256, 256) pixels, which yielded in this case an effective real
space pixel size of 15.3 nm. The dataset was acquired in step-scan mode with a
nominal scan window of 1 $\mu$m x 1 $\mu$m with a step size of 50 nm. This
gives 21 scan positions in vertical and horizontal scan direction, so that the
resulting dataset consisted of 441 diffraction patterns in total, each with an
exposure time of 0.5 s. The approximate scan positions were obtained from the
encoders of the piezo scanner and were used as a direct input for
ptychographic phase retrieval. The position errors of the setup were not
considered here. In practice, the beam, specimen and encoders can shift
relative to each other.

The following processing schemes were implemented using this test dataset:

  * •

1000 iterations on the full dataset for reference (Figure 5).

  * •

Following the original line scan pattern in y direction, add 10 patterns and
perform 10 iterations per step. After adding the complete data, iterate on the
full dataset until 1000 iterations in total are reached (Figure 6).

  * •

Add a random pick of 10 patterns and perform 10 iterations per step. Iterate
on the full dataset until 1000 iterations in total are reached (Figure 7).

  * •

Add 10 patterns following an inwards square spirangle pattern and perform 10
iterations per step. Iterate on the full dataset until 1000 iterations in
total are reached (Figure 8).

  * •

Split the dataset into four random disjoint subsets. Perform 10 iterations on
each subset independently, synchronize and repeat until 1000 iterations are
reached. This simulates parallel processing on offline data (Figure 9).

  * •

Add 10 patterns following the original line scan in y direction to four
disjoint subsets in a round-robin (Kleinrock, 1964) fashion, perform 10
iterations on each subset after adding patterns to a subset, and synchronize.
Iterate until 1000 total iterations are reached. This simulates parallel
processing on live data (Figure 10).

  * •

Following the original line scan pattern in y direction, add 10 patterns and
perform 1000 iterations per step to investigate impact of early over-
optimization (Figure 11).

See section “Data and code availability” for the test data, ePIE
implementation and test code.

###  2.6 Live ePIE

In the simulated live processing described above, the processing follows a
predefined sequence of data loading, iteration and visualization that is
executed sequentially in a deterministic fashion. All necessary parameters are
read from configuration files.

Figure 3: Diagram of the data flow for live ePIE.

After feasibility of live processing was demonstrated, these processes were
implemented for real-world live processing at the DESY beamline in a Jupyter
notebook as concurrent threads that follow the speed of the data source and
iterate as fast as possible (Figure 3). It was executed on a node of the
Maxwell cluster at DESY which has access to a cluster file system where
results from the beamline are written in chunks with a short time delay, as
well as the beamline’s event interface. In principle, the data could also be
streamed over network instead of using chunks in a file system, but such an
interface is currently not implemented at the beamline. The pipeline is armed
prior to data acquisition and starts processing as soon as a scan_started
event is received. Alternatively, the processing can attach to a running scan
or load data from previous scans. The expected scan area is extracted from the
scan_started event, while other reconstruction parameters are determined from
a configuration file for the ePIE engine.

A custom LiberTEM acquisition object was developed that reads the detector
data from the folder structure that the beamline creates. This acquisition
object is used as a data source for the preprocessing UDF. In addition,
another UDF reads the matching stage encoder values from the beamline data to
determine the position of each scan point. LiberTEM runs these UDFs in a pool
of worker processes in the background. The raw beamline data is available in
(Weber et al., 2023).

The run_udf_iter interface of LiberTEM provides a Python generator with
updated partial preprocessing result from the UDF for each received input data
partition, together with a map that describes which parts of the expected data
have been received so far. A feeder thread takes these updates from the
generator and feeds them into the processing thread.

The processing thread consumes these updates and waits until a sufficient
subset of data is received. Then it loads the data into the ePIE
implementation and starts iterating. Between iterations it checks if new data
is available and loads it as required, and regularly feeds an intermediate
state of the reconstruction into a queue for display.

The visualization thread reads from this display queue and updates LiberTEM
live plots with each new result.

Live processing was performed while recording data on an NTT-AT resolution
test chart with 50nm thinnest lines and spaces (model ATN/XRESO-50HC) at an
X-ray energy of E = 9 keV. The X-ray beam was focused by a Fresnel zone plate
(FZP) with outer diameter of 125 $\mu$m and thinnest outermost lines of 70 nm.
At the employed X-ray energy, these optics create a focus with a size of about
85 nm at a focal distance of 64 mm. The diffraction patterns were collected
with an Eiger X 4M detector (Dectris AG) positioned at a distance of 3.265 m
behind the sample. In the reconstructions, these patterns were cropped to a
size of (256, 256) pixels. The dataset was acquired using a continuous meander
scan with a nominal scan window of 10 $\mu$m x 11 $\mu$m and a step size of
100 nm. The resulting dataset consisted of 11245 diffraction patterns in
total, each with an exposure time of 20 ms.

The test object was mounted in a different orientation between the measurement
campaign for the simulated live processing and the actual live processing.
Both the recording for simulated live processing and for real live processing
were not adapted specifically for this purpose, but performed alongside
regular beamtimes at PtyNAMi during the preparation phase for real
measurements. The parameters are different and in the typical range for
PtyNAMi for that reason.

##  3 Results

###  3.1 Simulated live and parallel processing

Figure 4: Starting value for all simulated live reconstructions. The object is
initialized with the identity function. The illumination amplitude is
initialized with a small Gaussian peak and flat phase of 0. Figure 5:
Reference result with 1000 iterations on the full dataset. The magenta dots in
the upper left plot indicate the scan positions.

Figure 4 shows the starting value of illumination and object used in all
simulated live processing.

Figure 5 shows amplitude and phase of the reconstruction of both illumination
and object using 1000 iterations on the complete dataset, corresponding to
conventional offline processing.

Figures 6 to 10 show how the simulated live and parallel processing results
evolve with different processing and subdivision schemes. Note how the
reconstruction extends far beyond the scan range due to the extended size of
the illumination.

In all cases a reconstruction with strong resemblance to the final result
emerges already after a few iterations on a small subset of the dataset. The
random (Figure 7) and spirangle (Figure 8) scan patterns cover the entire scan
area early on, demonstrating how such acquisition schemes can be advantageous
for live processing compared to a sequential line scan to get an overview of
the object early on.

After completing the same number of iterations, the reconstructions that are
not parallelized differ from the reference result mainly by a small shift that
manifests itself as a contour line around sharp features of the object in the
difference plot. The parallelized results exhibit a stronger difference that
did not vanish up to 1000 iterations. The plot of amplitude and phase
difference in Figure 12 reveals that the amplitude of the object is larger and
the amplitude of the illumination is lower than the reference. However, they
still show a qualitative visual resemblance to the reference result.

Figure 6: Simulated live reconstruction following the original line scan
pattern in y direction. The magenta dots in the left column indicate the scan
positions. The right column shows the absolute value of the difference to the
reference result. Figure 7: Simulated live reconstruction with a random scan
pattern. The magenta dots in the left column indicate the scan positions. The
right column shows the absolute value of the difference to the reference
result. Figure 8: Simulated live reconstruction following a spirangle scan
pattern. The magenta dots in the left column indicate the scan positions. The
right column shows the absolute value of the difference to the reference
result. Figure 9: Simulated parallelized reconstruction on four random
disjoint subsets. The colored dots in the left column indicate the scan
positions, with the color indicating the subset. The right column shows the
absolute value of the difference to the reference result. Figure 10: Simulated
parallelized live reconstruction on four gradually growing disjoint subsets.
The colored dots in the left column indicate the scan positions, with the
color indicating the subset. The right column shows the absolute value of the
difference to the reference result. Figure 11: Simulated live reconstruction
following the original line scan pattern in y direction with over-optimization
on a small subset of input data. The magenta dots in the left column indicate
the scan positions. The right column shows the absolute value of the
difference to the reference result. Figure 12: Amplitude and phase difference
between the final result of Figure 10 and the reference result (Figure 5). The
amplitude of the object is higher, while the amplitude of the illumination is
lower. The contour line around the sharp features of the object phase indicate
a small shift relative to the reference result.

Animated videos of the simulated live reconstruction are available in (Weber
et al., 2023).

###  3.2 Live processing

Figure 13 shows how the live view of object and illumination reconstruction
develops during live processing. This allows monitoring the scan. Since the
reconstruction can be restarted and re-attached to a running scan, users can
already start optimizing reconstruction parameters while a scan is running.
See “Data and code availability” regarding availability of source code, raw
data and animated results.

Figure 13: Evolution of amplitude and phase for object and illumination that
was recorded during live processing. The screen with the Jupyter notebook
running the reconstruction was recorded as a video. This Figure shows
rearranged and annotated snapshots of the plots taken at the indicated times
relative to the start of the scan. The full video is available in the
Supplementary Material. The insets show the center of the illumination scaled
by a factor of 2. The first step shows the initialization before iteration
starts. The second step was recorded just after the start of iteration. Steps
3 and 4 show arbitrary intermediate states. The second to last step shows the
state at the moment the scan completes. The last snapshot was taken after a
number of additional iterations. In the beginning, 50 iterations were executed
for each batch of 500 images received from the beamline. Towards the end,
about 6 iterations were executed per batch due to the much larger number of
scan points included in each iteration. The live plots were updated every two
iterations for quick visual feedback.

A screen capture video of the live reconstruction is available in the
supplementary material.

##  4 Discussion

The results on simulated and real-world live processing show that iterative
ptychography algorithms can be used for live reconstruction by starting the
iteration on a small subset of the data that grows as the data acquisition
proceeds. The state of object and illumination are carried forward. For this
application a ptychography engine should offer the following features:

  * •

An API that allows integration in a complex workflow, as opposed to stand-
alone programs.

  * •

Flexible and fast import and export of object, illumination and data via the
API.

  * •

Stepping the iteration to extract intermediate values or stop as required.

  * •

Options to encourage stable convergence, such as good starting values for
illumination and object, restrictions on value ranges and illumination shapes,
countering impact of noise and/or detector gain characteristics, masking or
cropping invalid pixels, etc.

  * •

Configuration of a stable, pre-determined size and shape of the object and
illumination as well as a relation between physical scan position and pixel
position in the reconstructed object based on the expected acquisition
parameters.

  * •

Separation of preprocessing from the iteration engine.

In particular when the beamline at DESY acquires data at a low rate, around 10
frames per second, it is very advantageous to obtain a view of the object
early on during an acquisition to abort and adjust parameters in case issues
are detected, such as the choice of acquisition area, exposure, or focal
plane. During real-world live processing the implementation could easily keep
up with this data rate. During the experiment shown in Figure 13, the region
on the test specimen was already clearly recognizable after completion of
about a quarter of the scan. An operator has a chance to detect issues with
the field of view or settings at this point. Since large scans at PtyNAMi can
run for many hours overnight, detecting issues at 25% completion or earlier
can save significant amounts of beam time and allows early optimization of
reconstruction parameters.

Figure 11 shows that the reconstruction doesn’t diverge even with a high
number of iterations on a small initial subset for this particular dataset,
engine and choice of parameters. As soon as a region of the test specimen is
covered by a sufficient number of scan points, the iteration stabilizes on a
value similar to the reference result in this region. Similar to the other
scans, this reconstruction exhibits a small shift that results in a bright
line around sharp features in the difference plot. A result that clearly shows
features of the test object was obtained with three and a half scan rows in
this test. This indicates that a successful reconstruction requires some
degree of lateral coverage with overlapping illuminations, but starting the
iteration early with a smaller number of diffraction patterns is not harmful
in this case. The lower limit for successful partial reconstruction will
likely be different for different illumination shapes and sizes, scan
patterns, objects, and reconstruction engine settings. Typically, engines for
iterative ptychographic reconstruction will have many parameters and
constraints for regularization in order to ensure a stable convergence with
real-world data. Here, we merely show that the results for live processing for
all tests that we performed eventually converged towards the offline result as
more and more data was added.

As long as the convergence is stable, meaning it converges quickly towards the
same result independent of the starting values for illumination and object,
and stabilizes after a sufficient number of iterations, live reconstruction
will eventually converge on the offline result under all circumstances: The
previous reconstruction steps with a subset of the input data merely generate
a starting value for illumination and object for the final iteration phase on
the complete input data. The final iteration phase will then converge and
stabilize on the offline result.

Furthermore, any real-world dataset for ptychographic reconstruction can be
understood as a partial dataset: The actual scan size and the size of the
object are only limited by practical concerns and can be, in principle,
arbitrarily large. That means any real-world dataset can be extended by just
continuing the scan beyond the initial borders, making the initial “full”
dataset a partial dataset. The reconstruction on a subset of a large dataset
is not distinguishable from the reconstruction of a “full” dataset with the
same scan coordinates as the subset. Regions of the object transfer function
that are not touched by the illumination function for any scan position just
stay at their initial values.

The results shown here conform to this mathematical expectation: Eventually, a
reconstruction on a growing subset of a sufficiently large dataset will give a
stable reconstruction for all regions of the object that are covered by the
illumination function in that subset, provided the reconstruction of the
complete dataset is stable. For a general discussion and further references on
the convergence of ptychographic reconstruction see, for example, (Melnyk,
2023).

Simulated live reconstruction and actual live reconstruction mainly differ in
their input and output method: In simulated live processing, a “recipe” with
well-defined steps is executed, and all intermediate results are stored and
then plotted. In actual live processing, the sequence of data and iteration is
determined by the interplay of data source, preprocessing and iteration
thread. The intermediate results are not recorded here, but just plotted for
quick visual feedback. Notably, here live reconstruction is not meant to
replace careful offline reconstruction with optimized parameters. It is rather
a quick preview tool to observe what data is being recorded and to check if
the chosen parameters are likely to give a successful offline reconstruction
later.

In this paper, we show practical evidence that such live reconstruction can be
possible and beneficial under the typical conditions at PtyNAMi. In
particular, live reconstruction worked reliably for different illumination
shapes and scan parameters that were chosen somewhat arbitrarily from typical
values. The live reconstruction seems rather robust against variations in the
set of diffraction patterns and number of iterations: By adding more and more
diffraction patterns, the result reliably converged to closely resemble the
offline reconstruction result in the region covered by the scan, meaning it
seems to meet the conditions for successful partial reconstruction described
in the previous paragraphs. It seems intuitive that a scan that yields a
satisfactory reconstruction with only a subset of scan points is likely to
also yield a satisfactory result for the full scan. That means it is valuable
as a preview tool to give an early indication to operators if the scan and
reconstruction parameters will give a good final result. Beamline operator
feedback was positive.

The parallelized reconstruction schemes can potentially be used to generate
faster previews by using multiple GPUs in parallel. An engine and processing
scheme that produces the same result independent of the number of parallel
GPUs is required, which was not yet achieved here. Generally, algorithms that
act in a parallel manner, such as relaxed averaged alternating reflections
(RAAR) or difference map (DM), are advantageous for this (Enders & Thibault,
2016). The PtyPy software package (Enders & Thibault, 2016) already offers
parallelized reconstruction and can potentially be used as a reconstruction
engine for live processing in a similar fashion.

##  5 Outlook

Faster processing seems feasible, in particular targeting electron
ptychography with fast direct detection detectors and fourth generation
synchrotron radiation sources such as PETRA IV, since the preprocessing with
LiberTEM is efficient and highly scalable, and parallelized live processing on
several GPUs seems feasible with a suitable reconstruction engine. Further
optimization of the ePIE engine is thinkable as well, for example to run on
different platforms such as field-programmable gate arrays (FPGAs) using
Alpaka (Matthes et al., 2017).

In the future, live reconstruction can be factored in when choosing scan
patterns. A spiral or spirangle seems promising since it combines short
movements of the stage between subsequent acquisitions with early coverage of
the whole field of view, and such patterns are already available at PtyNAMi. A
conventional line scan pattern can be advantageous if the position encoders
exhibit hysteresis. Further experiments to test the behavior of live
reconstruction for different scan patterns, specimens, illuminations and scan
step sizes seem promising. Here, we demonstrate that live reconstruction can
yield useful results for different scan patterns.

The approach for live iterative ptychography described here is not limited to
ePIE, but should work with any reconstruction engine that offers a similar
interface. Implementations of different reconstruction approaches such as
multislice ptychography, mixed state ptychography, near field ptychography
(You et al., 2023) and/or alternative projection schemes such as RAAR,
alternating direction method of multipliers (ADMM) and similar can meet the
requirements for live processing described in the previous section: They
optimize gradually, and the set of scan positions is flexible. Furthermore,
this approach should be transferable to other experimental modalities such as
extreme ultraviolet light, visible light and electrons. A trial with electron
microscopy data is included in the Supplementary Material.

Live monitoring is not limited to iterative ptychographic reconstruction: Any
other processing or reconstruction method that is compatible with LiberTEM can
be included with little effort in the preprocessing stage, such as Centre of
Mass or even direct ptychographic reconstruction (Strauch et al., 2021; Bangun
et al., 2022).

Furthermore, in the future it can be investigated under which conditions the
reconstruction is quantitative and how large the errors are at various stages
of the scan and reconstruction.

Panning and zooming with a live view may be possible by applying a suitable
transformation to an intermediate state of the object in order to generate a
new starting value, generating a new configuration and continuing the
acquisition and iteration with the new parameters. This can allow live
navigation on the specimen and adjustment of parameters at the beam line,
effectively creating an interactive microscope for X-rays and other
ptychography modalities. The results also show that non-standard scan patterns
can optimize early coverage of the field of view. Previously, Velazco et al.
(Velazco et al., 2020) have shown non-standard scan strategies, including a
space-filling curve. Live reconstruction can also enable adaptive scanning
strategies (Ede, 2021).

ePIE

    extended ptychographical iterative engine
REPL

    read–eval–print loop
UDF

    user-defined function
SSB

    single side-band
WDD

    Wigner distribution deconvolution
ADMM

    alternating direction method of multipliers
RAAR

    relaxed averaged alternating reflections
DM

    difference map
DFT

    Discrete Fourier Transform

##  6 Data and code availability

The ePIE reconstruction engine and interface for live processing is available
upon reasonable request from Andreas Schropp<EMAIL_ADDRESS>
Additional movies that show the live reconstruction, the dataset and notebook
for simulated offline processing, the Jupyter notebook for live processing at
the beamline, and the beamline data recorded during the live processing run
are available at http://doi.org/10.5281/ZENODO.8239052 (Weber et al., 2023).

##  7 Competing interests

The authors declare no competing interests.

##  8 Acknowledgements

  * •

Parts of this research were carried out at beamline P06 of PETRA III at
Deutsches Elektronen-Synchrotron DESY, a member of the Helmholtz Association
(HGF).

  * •

We acknowledge support from Helmholtz Association under contract No. ZT-I-0025
(Ptychography 4.0) and JL-MDMC (Joint Lab on Model and Data-Driven Material
Characterization).

  * •

This research was supported in part through the Maxwell computational
resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany

  * •

This work was partly funded by the Center for Advanced Systems Understanding
(CASUS) which is financed by Germany’s Federal Ministry of Education and
Research (BMBF) and by the Saxon Ministry for Science, Culture, and Tourism
(SMWK) with tax funds on the basis of the budget approved by the Saxon State
Parliament.

## References

  * Bangun et al. (2022) Bangun, A., Baumeister, P.F., Clausen, A., Weber, D. & Dunin-Borkowski, R.E. (2022). Wigner distribution deconvolution adaptation for live ptychography reconstruction, URL https://arxiv.org/abs/2212.01309.
  * Ede (2021) Ede, J.M. (2021). Adaptive partial scanning transmission electron microscopy with reinforcement learning, Machine Learning: Science and Technology 2, 045011.
  * Enders & Thibault (2016) Enders, B. & Thibault, P. (2016). A computational framework for ptychographic reconstructions, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472, 20160640.
  * Hoppe (1969a) Hoppe, W. (1969a). Beugung im inhomogenen Primärstrahlwellenfeld. I. Prinzip einer Phasenmessung von Elektronenbeungungsinterferenzen, Acta Crystallographica Section A 25, 495–501, URL https://doi.org/10.1107/S0567739469001045.
  * Hoppe (1969b) Hoppe, W. (1969b). Beugung im inhomogenen Primärstrahlwellenfeld. III. Amplituden- und Phasenbestimmung bei unperiodischen Objekten, Acta Crystallographica Section A 25, 508–514, URL https://doi.org/10.1107/S0567739469001069.
  * Hoppe & Strube (1969) Hoppe, W. & Strube, G. (1969). Beugung in inhomogenen Primärstrahlenwellenfeld. II. Lichtoptische Analogieversuche zur Phasenmessung von Gitterinterferenzen, Acta Crystallographica Section A 25, 502–507, URL https://doi.org/10.1107/S0567739469001057.
  * Kleinrock (1964) Kleinrock, L. (1964). Analysis of a time‐shared processor, Naval Research Logistics Quarterly 11, 59–73.
  * Könnecke et al. (2015) Könnecke, M., Akeroyd, F.A., Bernstein, H.J., Brewster, A.S., Campbell, S.I., Clausen, B., Cottrell, S., Hoffmann, J.U., Jemian, P.R., Männicke, D., Osborn, R., Peterson, P.F., Richter, T., Suzuki, J., Watts, B., Wintersberger, E. & Wuttke, J. (2015). The nexus data format, Journal of Applied Crystallography 48, 301–305.
  * Maiden & Rodenburg (2009) Maiden, A.M. & Rodenburg, J.M. (2009). An improved ptychographical phase retrieval algorithm for diffractive imaging, Ultramicroscopy 109, 1256–1262, URL https://www.sciencedirect.com/science/article/pii/S0304399109001284.
  * Matthes et al. (2017) Matthes, A., Widera, R., Zenker, E., Worpitz, B., Huebl, A. & Bussmann, M. (2017). Tuning and optimization for a variety of many-core architectures without changing a single line of implementation code using the alpaka library, J.M. Kunkel, R. Yokota, M. Taufer & J. Shalf (eds.), High Performance Computing, 496–514, Cham: Springer International Publishing.
  * Melnyk (2023) Melnyk, O. (2023). Convergence properties of gradient methods for blind ptychography, arXiv .
  * Nashed et al. (2014) Nashed, Y.S.G., Vine, D.J., Peterka, T., Deng, J., Ross, R. & Jacobsen, C. (2014). Parallel ptychographic reconstruction, Optics Express 22, 32082.
  * Pelz et al. (2021) Pelz, P.M., Johnson, I., Ophus, C., Ercius, P. & Scott, M.C. (2021). Real-time interactive 4D-STEM phase-contrast imaging from electron event representation data: Less computation with the right representation, IEEE Signal Processing Magazine 39, 25–31.
  * Rodenburg & Faulkner (2004) Rodenburg, J.M. & Faulkner, H.M. (2004). A phase retrieval algorithm for shifting illumination, Applied physics letters 85, 4795–4797.
  * Schroer et al. (2005) Schroer, C.G., Kurapova, O., Patommel, J., Boye, P., Feldkamp, J., Lengeler, B., Burghammer, M., Riekel, C., Vincze, L., van der Hart, A. & Kuchler, M. (2005). Hard x-ray nanoprobe based on refractive x-ray lenses, Applied Physics Letters 87, 124103, URL http://apl.aip.org/applab/v87/i12/p124103_s1.
  * Schropp et al. (2010) Schropp, A., Boye, P., Feldkamp, J.M., Hoppe, R., Patommel, J., Samberg, D., Stephan, S., Giewekemeyer, K., Wilke, R.N., Salditt, T., Gulden, J., Mancuso, A.P., Vartanyants, I.A., Weckert, E., Schöder, S., Burghammer, M. & Schroer, C.G. (2010). Hard x-ray nanobeam characterization by coherent diffraction microscopy, Applied Physics Letters 96, 091102, URL http://link.aip.org/link/?APPLAB/96/091102/1http://aip.scitation.org/doi/10.1063/1.3332591, nULL.
  * Schropp et al. (2020) Schropp, A., Döhrmann, R., Botta, S., Brückner, D., Kahnt, M., Lyubomirskiy, M., Ossig, C., Scholz, M., Seyrich, M., Stuckelberger, M.E., Wiljes, P., Wittwer, F., Garrevoet, J., Falkenberg, G., Fam, Y., Sheppard, T.L., Grunwaldt, J.D. & Schroer, C.G. (2020). PtyNAMi: ptychographic nano-analytical microscope, Journal of Applied Crystallography 53, 957–971, URL https://doi.org/10.1107/S1600576720008420.
  * Schropp & Schroer (2010) Schropp, A. & Schroer, C.G. (2010). Dose requirements for resolving a given feature in an object by coherent x-ray diffraction imaging, New Journal of Physics 12, 035016, URL https://dx.doi.org/10.1088/1367-2630/12/3/035016.
  * Strauch et al. (2021) Strauch, A., Weber, D., Clausen, A., Lesnichaia, A., Bangun, A., März, B., Lyu, F.J., Chen, Q., Rosenauer, A., Dunin-Borkowski, R. & et al. (2021). Live processing of momentum-resolved STEM data for first moment imaging and ptychography, Microscopy and Microanalysis 1–15.
  * Velazco et al. (2020) Velazco, A., Nord, M., Béché, A. & Verbeeck, J. (2020). Evaluation of different rectangular scan strategies for STEM imaging, Ultramicroscopy 215, 113021, URL https://www.sciencedirect.com/science/article/pii/S0304399120300334.
  * Weber et al. (2023) Weber, D., Ehrig, S., Schropp, A., Clausen, A., Achilles, S., Hoffmann, N., Bussmann, M., Dunin-Borkowski, R.E. & Schroer, C.G. (2023). Live iterative ptychography, URL http://doi.org/10.5281/ZENODO.8239052.
  * You et al. (2023) You, S., Lu, P.H., Schachinger, T., Kovács, A., Dunin-Borkowski, R.E. & Maiden, A.M. (2023). Lorentz near-field electron ptychography, Applied Physics Letters 123.

  *[SSB]: single side-band
  *[WDD]: Wigner distribution deconvolution
  *[ePIE]: extended ptychographical iterative engine
  *[DFT]: Discrete Fourier Transform
  *[REPL]: read–eval–print loop
  *[UDF]: user-defined function
  *[RAAR]: relaxed averaged alternating reflections
  *[DM]: difference map
  *[ADMM]: alternating direction method of multipliers


#  The dust sublimation region of the Type 1 AGN NGC4151 at a hundred micro-
arcsecond scale
as resolved by the CHARA Array interferometer

Makoto Kishimoto  Department of Astrophysics & Atmospheric Sciences, Kyoto
Sangyo University, Kamigamo-motoyama, Kita-ku, Kyoto 603-8555, Japan  Matt
Anderson  Georgia State University Center for High Angular Resolution
Astronomy Array, Mount Wilson Observatory, USA  Theo ten Brummelaar  Georgia
State University Center for High Angular Resolution Astronomy Array, Mount
Wilson Observatory, USA  Christopher Farrington  Georgia State University
Center for High Angular Resolution Astronomy Array, Mount Wilson Observatory,
USA  Robert Antonucci  Department of Physics, University of California, Santa
Barbara, USA  Sebastian Hoenig  Department of Physics & Astronomy, University
of Southampton, UK  Florentin Millour  Laboratoire Lagrange, Université Côte
d’Azur, Observatoire de la Côte d’Azur, CNRS, Nice, France  Konrad Tristram
European Southern Observatory, Santiago, Chile  Gerd Weigelt  Max Planck
Institute for Radio Astronomy, Bonn, Germany  Laszlo Sturmann  Georgia State
University Center for High Angular Resolution Astronomy Array, Mount Wilson
Observatory, USA  Judit Sturmann  Georgia State University Center for High
Angular Resolution Astronomy Array, Mount Wilson Observatory, USA  Gail
Schaefer  Georgia State University Center for High Angular Resolution
Astronomy Array, Mount Wilson Observatory, USA  Nic Scott  Georgia State
University Center for High Angular Resolution Astronomy Array, Mount Wilson
Observatory, USA

###### Abstract

The nuclear region of Type 1 AGNs has only been partially resolved so far in
the near-infrared (IR) where we expect to see the dust sublimation region and
the nucleus directly without obscuration. Here we present the near-IR
interferometric observation of the brightest Type 1 AGN NGC4151 at long
baselines of $\sim$250 m using the CHARA Array, reaching structures at hundred
micro-arcsecond scales. The squared visibilities decrease down to as low as
$\sim$0.25, definitely showing that the structure is resolved. Furthermore,
combining with the previous visibility measurements at shorter baselines but
at different position angles, we show that the structure is elongated
perpendicular to the polar axis of the nucleus, as defined by optical
polarization and a linear radio jet. A thin-ring fit gives a minor/major axis
ratio of $\sim$0.7 at a radius $\sim$0.5 mas ($\sim$0.03 pc). This is
consistent with the case where the sublimating dust grains are distributed
preferentially in an equatorial plane in a ring-like geometry, viewed at an
inclination angle of $\sim$40°. Recent mid-IR interferometric finding of
polar-elongated geometry at a pc scale, together with a larger-scale polar
outflow as spectrally resolved by the HST, would generally suggest a dusty,
conical and hollow outflow being launched presumably in the dust sublimation
region. This might potentially lead to a polar-elongated morphology in the
near-IR, as opposed to the results here. We discuss a possible scenario where
an episodic, one-off anisotropic acceleration formed a polar-fast and
equatorially-slow velocity distribution, having lead to an effectively flaring
geometry as we observe.

Active galactic nuclei — Interferometry

##  1 Introduction

With long-baseline infrared interferometry, it is possible to spatially
resolve the thermal emission from nuclear structures on sub-parsec scales in
the nuclei of active galaxies (Active Galactic Nuclei, or AGN). Such AGN are
thought to comprise a supermassive black hole surrounded by optically thick
material that obscures sight of the central putative accretion disk and
ionized gas if viewed from equatorial direction, while this nucleus is
directly seen along polar direction (Antonucci & Miller, 1985). The latter
polar-viewed objects are believed to be Type 1s, while the former
equatorially-viewed ones are thought to be Type 2s. The system of this
optically-thick, dusty material has been called a ’torus’.

Mid-infrared (mid-IR) interferometry explores the relatively outer region of
the torus where the temperature of dust grains heated by the nucleus is
$\sim$300 K, while the interferometry in the near-IR can probe the innermost
dusty structure where dust grains are almost at the sublimation temperature of
$\sim$1500 K. Here in this paper, we focus on this innermost region. For its
exploration, the most suitable objects would be Type 1 AGNs, where the
innermost region is directly seen without significant obscuration and without
much complication and confusion from dust extinction.

However, the overall angular size of the Type 1 nuclei turned out to be quite
small. The innermost dusty region, or the dust sublimation region, is expected
to be larger for brighter objects (the angular size would roughly be
proportional to the square-root of apparent flux; e.g. Figure 5 of Kishimoto
et al. 2007), but even in the brightest Type 1 source on the sky, NGC4151, the
region is only partially resolved with the Keck interferometer’s (hereafter
KI) 85 m baseline (i.e. shows only a small decrease of visibility; Swain et
al. 2003; Kishimoto et al. 2009; Pott et al. 2010). Other Type 1 AGNs are also
only partially resolved with Very Large Telescope Interferometer (VLTI)
baselines up to 130 m (Weigelt et al. 2012; GRAVITY Collaboration et al. 2019;
Leftley et al. 2021).

Another way to measure the radius of the dust sublimation region is to detect
the time lag of the near-IR emission variation with respect to that of the
nuclear optical emission (Koshida et al. 2009 and references therein), which
has been very successful. The interferometric radii inferred from the small
visibility decrease described above are in fact consistent with the time-lag
radius where the interferometric radii, which are brightness-weighted radii
giving the overall size of the dust sublimation region, are slightly larger
than the near-IR reverberation radii, which are response-weighted providing
the inner boundary radius (e.g. Kishimoto et al. 2009, 2011a; for NGC4151, the
interferometric radius is $\sim$0.5 mas, while the time-lag radius corresponds
to $\sim$0.3 mas; see Figures 3 and 4 in Kishimoto et al. 2013). If this
interferometric size estimate is really correct, the observations with much
longer baselines would reveal a large decrease of visibility, showing a
definitely resolved structure. Here we present observations of NGC4151 with
the CHARA Array at $\sim$250 m baselines, a factor of $\sim$2 longer than
previously achieved for AGNs, or any other extragalactic objects in fact. We
show that the near-IR visibility goes down significantly at this long baseline
and argue that we started to see morphological information on this target.

##  2 Observations

We have observed NGC4151 with the CHARA Array (ten Brummelaar et al., 2005)
over several different nights as summarized in Table 1. We used two baselines
– a short, $\sim$34 m baseline using the telescopes S1 and S2 in 2020, and a
long, $\sim$250 m baseline with telescopes S2 and W1 in 2021. In all
occasions, the telescopes were equipped with an adaptive optics (AO) system,
which has been commissioned recently. The details and operations of the CHARA
Array AO system are described by ten Brummelaar et al. (2018) and Anugu et al.
(2020), while the early design of the AO system is described by Che et al.
(2013). We used the instrument CLASSIC, a two-beam combiner which is the most
sensitive to faint targets among the various instruments at the CHARA Array.
We used three different calibrators listed in Table 2, which are very close in
projected distance to NGC4151. They vary in optical and infrared magnitudes,
with TIC9453648 closest in optical magnitude to NGC4151. The science target
was bracketed by calibrator observations close in time (except for the
observation on 2020-02-15 when the time interval for one of the calibrator
observations was long; see Appendix). We observed them all in K’-band
($\lambda_{\rm obs}$ = 2.13 $\mu$m), not spectrally dispersed. The details of
the sequences and raw visibilities observed are shown in Appendix.

Table 1: Observed visibilities for NGC4151 at CHARA/CLASSIC in K’ band Observing | Time | telescopes | projected | baseline | $V^{2}$ |  $V^{2}$ corr. for
---|---|---|---|---|---|---
Date (UT) | (UT) |  | baseline (m) | PA (°) |  |  $f_{\rm AD}$=0.16$\pm$0.05
2020-02-15 | 09h21m | S1-S2 | 33.9 | 179.0 | 0.87 $\pm$ 0.15 | 0.85 $\pm$ 0.17
2020-02-16 | 09h21m | S1-S2 | 33.9 | 7.3 | 1.00 $\pm$ 0.27 | 1.00 $\pm$ 0.33
2021-03-19 | 07h44m | S2-W1 | 249.1 | 142.2 | 0.25 $\pm$ 0.08 | 0.17 $\pm$ 0.08
2021-04-30 | 05h18m | S2-W1 | 249.4 | 139.1 | 0.36 $\pm$ 0.09 | 0.28 $\pm$ 0.10
2021-04-30 | 08h30m | S2-W1 | 216.6 | 113.3 | 0.63 $\pm$ 0.07 | 0.57 $\pm$ 0.08
2021-04-30 | 09h23m | S2-W1 | 193.2 | 107.4 | 0.57 $\pm$ 0.15 | 0.50 $\pm$ 0.16
Table 2: Calibrators and projected distance $r_{\rm proj}$ from NGC4151.

name | V (mag) | K (mag) |  $r_{\rm proj}$ (°) | diameter (mas)a
---|---|---|---|---
HD105881 | 8.0 | 7.0 | 0.3 | 0.161 $\pm$ 0.004
SAO62878 | 9.8 | 8.3 | 0.09 | 0.092 $\pm$ 0.002
TIC9453648 | 11.6 | 8.7 | 0.2 | 0.090 $\pm$ 0.002

a Diameters from SearchCal by JMMC.

##  3 Data processing

The data were processed following the steps described by Ten Brummelaar et al.
(2013). The instrument CLASSIC scans each fringe with a dithering mirror. For
NGC4151 and all calibrators, we needed to scan relatively slowly, taking about
1.5 seconds to scan over $\sim$70 $\lambda_{\rm obs}$ (scanning frequency of
$\sim$250 Hz, 5 samples per fringe, scanning over 150 $\mu$m). After
subtracting an averaged sky scan (scan of nearby sky region which is regarded
as a ‘dark’ frame) from each scan, we first normalize each target scan by its
low-pass-filtered version. Then we Fourier-transform the normalized scan to
measure its fringe power by computing the power spectrum. The noise bias has
been subtracted using the average power spectrum of the off-fringe scans.

We have constructed our own IDL scripts to implement these steps, and
calculated 10-scan, or 15-second, averages of squared visibilities $V^{2}$
over each fringe track which lasted about 10 minutes (5 minutes on fringe + 5
minutes off fringe). This confirmed overall stability of the fringes, but also
clearly identified certain tracks where fringes were within the scan range
only in some portion of the track. In these cases, we only used one continuous
range per track where fringes are visible in the $V^{2}$ track. Then we
integrated all the power spectra from each track (or the designated range of
the track) to calculate the final raw $V^{2}$ and also its uncertainty from
the fluctuation of the raw $V^{2}$ over the track. The fluctuation of the off-
fringe scans is also calculated and added in quadrature to the raw $V^{2}$
uncertainty. The system visibility for each observation of NGC4151 was
determined from the adjacent measurements of calibrators. Finally we obtained
the final calibrated visibility by dividing the raw $V^{2}$ with the estimated
system visibility. The details of these processing and raw $V^{2}$ tracks are
shown in Appendix.

Figure 1: Observed visibility squared $V^{2}$ as a function of projected
baselines. In addition to our CHARA Array data at $\sim$30 m and
$\sim$200$-$250 m baselines, we also plotted visibilities observed at Keck
interferometer (Kishimoto et al. 2013 and references therein). $V^{2}$ for
gaussian with HWHM of 0.3 and 0.4 mas are also plotted for reference in dotted
and dashed curves, respectively. Figure 2: Sampled uv points, color-coded with
observed visibility squared.

##  4 Results

###  4.1 Visibility vs baseline

Figure 1 shows all the calibrated squared visibilities $V^{2}$ measured at
short ($\sim$34 m) and long projected baselines ($\sim$190$-$250 m) with the
CHARA Array, together with the previous measurements using the Keck
interferometer (KI; Swain et al. 2003; Kishimoto et al. 2009; Pott et al.
2010; Kishimoto et al. 2011a, 2013). 111We have VLTI/AMBER data for NGC4151
taken in 2012 as shown in Kishimoto et al. (2013), but we do not include them
in our analysis here, since the VLTI/AMBER data might be affected by the
coherence loss due to AO performance which we now fear has not been well
quantified. However, the radius measurement for 2012 is represented by our KI
measurements in 2012, which we include here. The radius implied by the
combination of KI and VLTI data is consistent with the one by KI, as
quantified in Table 1 of Kishimoto et al. (2013). The KI measurements only
showed a hint of $V^{2}$ decrease as a function of baselines with only a small
range of projected baseline lengths. Thus the data did not rule out the case
where there is an unresolved source of a high flux fraction ($\sim$90%) with
the rest of the flux coming from a large scale structure. On the other hand,
for different Type 1 AGNs, VLTI instruments AMBER and GRAVITY have shown that
the visibilities decrease slightly over a wider range of baselines from
$\sim$40 m to 130 m (Weigelt et al. 2012; GRAVITY Collaboration et al. 2019).
Here we show that the visibility does actually decrease significantly at long
baselines in NGC4151, establishing that we are definitely resolving the
structure spatially.

###  4.2 Position angle dependence

The CHARA Array measurements reported here and the existing KI measurements
for NGC4151 are in fact complementary in the coverage of the position angle
(PA) of the projected baselines. Figure 2 shows the uv coverage of all the
observations, color-coded with the measured $V^{2}$.

Firstly, the four data points at 190$-$250 m baselines obtained with the CHARA
Array actually suggest a dependency of $V^{2}$ on the PA, in the sense that
the structure seems to be a bit more compact at PA$\sim$110° than at
$\sim$140°. To further investigate this, we show in Figure 3 the implied thin-
ring radius for each $V^{2}$ measurement along the PA of the baselines for the
CHARA Array measurements, as well as the KI measurements (see more below for
the corrections implemented). The three $V^{2}$ data points which are
consistent with unity within uncertainties are not shown in this plot, since
they do not give meaningful radius constraints alone. As a whole, we seem to
see a moderate elongation along roughly a north-south direction.

Here we have chosen to use a thin-ring model rather than a Gaussian geometry,
since the central part is considered devoid of dust emission due to dust
sublimation (see more discussions in section 5.4). We are still unable to
differentiate these geometries from the data themselves which are within the
first lobe of the visibility curve. The ring is adopted to be infinitely thin
to minimize the number of parameters involved.

Furthermore, we do expect that the central part, even if devoid of dust
emission, has another thermal emission component from the putative accretion
disk (AD). Based on the optical to near-IR spectral energy distribution, this
has a fractional flux contribution at K-band $f_{\rm AD}$ of $\sim$$0.16\pm
0.05$ for NGC4151 (Kishimoto et al. 2013 and references therein). Therefore,
assuming that this component remains unresolved, we calculated the $V^{2}$ for
the dust emission component (quoted in Table 1) and the corresponding thin-
ring radii. These are the radii shown in Figure 3.

In Table 3, we show the results of the thin-ring fit to the whole $V^{2}$ data
set, both for the simple ring and AD-corrected ring fit (note the difference
is only $\sim$10%). The latter is shown in Figure 4$a$ in dotted line, on top
of the same ring radii as shown in Figure 3. The projected ring has its major
axis at a PA of 17$\pm$10° with a radius 0.507$\pm$0.047 mas. The minor to
major axis ratio is 0.74$\pm$0.09, thus showing an elongation.

The lack of baseline length coverage between 90 and 190 m leaves some
uncertainty in the geometric fit to reproduce the data. An alternative would
be a power-law brightness distribution with a certain inner cut-off, a
distribution that is more extended in the outer region and more compact in the
inner region. However, we found that such a shallow power-law fit to the
observed visibilities yields a too small inner cut-off radius, which would be
inconsistent with the measured time-lag radius (see section 1). Therefore
despite the baseline gap, the adopted thin-ring model seems quite adequate in
measuring the overall size and elongation direction of the brightness
distribution.

Figure 3: Thin-ring radius implied by each visibility observed (after the
small correction for the putative AD component), plotted along the PA of the
corresponding baseline. The PA of the optical polarization (94°), indicative
of the system polar axis, is shown as the inner dashed-solid line, while the
PA of the linear radio jet structure (84°; components D1/D2/D3 over 50 mas
scale; Mundell et al. 2003) is indicated with the outer dashed-solid line
segments.  Table 3: Thin-ring fit results

| radius (mas) | min/maj | PA (°) |
---|---|---|---|---
all data | $0.462\pm 0.050$ | $0.74\pm 0.10$ | $19\pm 10$ |
CHARA + 2012 KI | $0.669\pm 0.086$ | $0.46\pm 0.08$ | $25\pm 6$ |
all data with | $0.507\pm 0.047$ | $0.74\pm 0.09$ | $17\pm 10$ |
$f_{\rm AD}$=0.16$\pm$0.05 |

For the KI data, we incorporate the conservative systematic uncertainty of
0.03 in $V^{2}$ (Section 3.2 in Kishimoto et al. 2011a).

Figure 4: ($a$) The same ring radii as in Figure 3 but with the fitted thin-
ring model shown in dotted line. The inner dashed-solid line indicates the
polar axis as probed by the optical polarization, and the outer dashed-solid
line segments show the PA of the radio jet, exactly the same as in Figure 3.
($b$) HST/WFPC2 archival image of the narrow-line region in [OIII]5007 line,
with the two sets of dashed-solid line segments corresponding to those in
panel $a$. ($c$) Three-dimensional distribution of [OIII]-emitting clouds
reconstructed from HST/STIS multi-slit data, which cover the volume indicated
by gray lines. The HST image from panel (b) in green is also shown at the back
of the 3D cube for reference, as the 2D flux distribution projected on to the
sky plane. This is as seen almost perpendicular to our line of sight. An
approximate axis of this 100-pc-scale conical structure is indicated in a
dashed-solid line.

##  5 Discussions

###  5.1 System axis direction and equatorial elongation

The system axis direction of the nucleus of NGC4151 can be inferred from a few
different observations. At a few arcsec scale, HST images show that the
narrow-line region is extended roughly along an E-W or NE-SW direction as
shown in Figure 4$b$, indicating that the polar axis of the system is along
this direction. A more accurate reference at a slightly smaller spatial scale
is the linear jet-like structure observed at radio wavelengths. The one-
arcsec-scale radio map clearly shows a linear structure roughly extended along
E-W direction, with some systematic wiggles along the structure (Mundell et
al., 2003). In the innermost region where the component D is considered to
have the active nucleus (Mundell et al. 2003; Ulvestad et al. 2005), the
structure is linearly extended along $\sim$84°over $\sim$50 mas scale at
$\sim$10 mas resolution (Mundell et al., 2003). The brightest sub-component D3
has been shown to be also extended roughly along the same direction at a
$\sim$1 mas resolution (Ulvestad et al. 2005; see section 5.6).

The reference for the axis direction at an even smaller scale is from optical
spectropolarimetry. NGC4151 shows an optical polarization of $\sim$1%, with
its PA parallel to the linear jet direction, which is quite typical for Type 1
AGNs (Antonucci 1983). The polarization is believed to originate from electron
scattering in an optically-thin, equatorial region (Antonucci 1983; Smith et
al. 2004; Lira et al. 2020), located somewhere between the broad-line region
(or co-spatial with the broad-line region) and dust sublimation region. Thus,
in terms of the spatial scale, this is probably the most suitable reference.
The polarization seems to be variable, namely 91° and 97°, observed $\sim$2
years apart (1992 / 1994-1995; Martel 1998)

These directions are graphically shown in Figure 3 – the inner dashed-solid
line indicates the polarization PA (at 94°, the average of the two above),
while the outer dashed-solid line segments corresponds to the linear jet
direction. The same directions are shown also in Figure 4$a$ and $b$. If we
consider a roughly axisymmetric structure, these directions are supposed to
indicate the polar axis direction of the system. Therefore, the elongation
observed with the near-IR interferometry (thin-ring fit is shown in Figure
4$a$) is approximately perpendicular to this system axis, implying that the
near-IR continuum emitting region looks elongated rather along the equatorial
direction of the system.

###  5.2 Size variability

Due to the variability of the nuclear luminosity, the radius of the dust
sublimation region could vary over the years. The radius of the near-IR
emitting region in NGC4151 has been claimed to show variation both in the
near-IR reverberation (Minezaki et al., 2004; Koshida et al., 2009; Schnülle
et al., 2015) and interferometry (Kishimoto et al., 2013). Pott et al. (2010)
pointed out no significant increase in interferometric radii against the flux
increase by a factor of 2, while Kishimoto et al. (2013) argued that the
radius change might be occurring as a function of the averaged AD flux over
several years. On the other hand, Schnülle et al. (2015) showed that the time-
lag radius decreased as opposed to the interferometric radius increase. If the
radius does vary, however, it might not be adequate to compare KI and CHARA
data here, since the KI and CHARA observations are apart in time over $\sim$10
years.

To quantify its potential effect, we have estimated the current expected size
based on the prescription for the variability given in Kishimoto et al.
(2013), where the size is not determined by the instantaneous nuclear flux,
but it varies as a function of its long-time average (because the size
variation would depend on the timescale of dust sublimation and dust
formation). We supplemented the light curve shown in Figure 3 of Kishimoto et
al. (2013) with the monitoring data available from
AAVSO222http://www.aavso.org for NGC4151 over the last $\sim$10 years, and
calculated the suggested 6-year flux average. The flux has been relatively
steady since the last high state occurred around 2011, and the resulting
predicted size at the current period is calculated to be roughly the same
(within $\lesssim 5$%) as that in 2012. We have then implemented a thin-ring
fit to the CHARA data only with the 2012 KI data (i.e. exclude the ones from
2009-2011). The results are shown in Table 3. The minor-to-major axis ratio
becomes smaller, but the elongation PA does not change significantly so that
the elongation direction still seems roughly perpendicular to the system axis.

###  5.3 Inclined ring and inclination angle

The near-IR continuum emission at $\sim$2 $\mu$m is believed to be
predominantly from the innermost, hottest dust grains almost at the
sublimation temperature. We tend to consider these sublimating dust grains to
be concentrated in the equatorial plane of the system, probably slightly
outside the broad-emission-line region. The elongation observed above is
consistent with such a naive expectation. In this case, we can formally derive
the inclination angle from the ellipticity of the elongation, assuming that
the region has a ring-like structure. From the minor to major axis ratio of
the thin-ring fit described above (Figure 4$a$), we would infer an inclination
angle of $\sim$43 $\pm$ 6°.

We can compare this value with other estimations from different data sets.
Among others, Das et al. (2005) deduced an inclination angle of 45$\pm$5° by
modeling the HST/STIS multi-slit data with a conical hollow outflow model over
a 100 pc scale region. The modeling is based on the assumption of a “linear”
outflow where the velocity of clouds are radial in direction and linearly
proportional to the distance from the nucleus. In fact, if we adopt the same
velocity field, the distribution of line-emitting material in three-
dimensional space can directly be reconstructed from the data, as shown by
Miyauchi & Kishimoto (2020) for a similar HST/STIS data set for the Seyfert 2
galaxy NGC1068. Here we have implemented the same 3D reconstruction using the
data set of Das et al. (2005) for NGC4151. The result is shown in Figure 4$c$.
Details are described elsewhere, but the 3D distribution visualizes and
illustrates the possible configuration of the material at a 100 pc scale, or 1
arcsec scale. Consistent with the modeling of Das et al. (2005), it shows a
patchy, but approximately hollow and conical distribution, and does suggest an
inclination consistent with the value inferred from the assumption of an
equatorial ring for the interferometry data. The SW side is the conical
outflow toward us, as deduced by Das et al. (2005). It is quite striking that
the 100 pc scale structure and the much smaller scale structure probed by our
CHARA Array observations suggest similar inclinations.

Figure 5:  ($a$) Schematic 3D geometry inferred for the central region
(reproduced from Figure 8b in Miyauchi & Kishimoto 2020). ($b$) 2D cut of the
3D geometry shown in panel $a$ (reproduced from Figure 8d in Miyauchi &
Kishimoto 2020 with additional annotations). ($c$) The same geometry as in
panel $a$, but rotated in 3D to match the possible configuration in NGC4151 as
seen from us, where the central dotted square schematically corresponds to the
region mapped by the CHARA Array and KI as indicated in panel $e$, which is
the same as Figure 3. ($d$) The same 3D geometry but further rotated to a side
view, perpendicular to our line of sight to the object. Approximate 3D
locations inferred for the VLBI components (15 GHz; Ulvestad et al. 2005) are
schematically shown in red circles. ($e$) See the description for panel $c$.
($f$) Reconstructed 3D distribution shown in Figure 4c, but with approximate
3D positions inferred for the knots A–H of the radio jet (Mundell et al.
2003). The molecular hydrogen distribution as inferred by Fernandez et al.
(1999) is shown as a red transparent ring in 3D with inner and outer radii of
0.5 and 1.5 arcsec.

###  5.4 Why in an equatorial plane?

In the mid-IR, long-baseline interferometry has shown that many objects have
polar-elongated morphology at pc scales, or tens of $R_{\rm sub}$ scales where
$R_{\rm sub}$ is the dust sublimation radius (Hönig et al. 2012, 2013;
Tristram et al. 2014; López-Gonzaga et al. 2014, 2016), though the mid-IR
interferometry data for NGC4151 are not clear about this point, since the PA
coverage of the data is poor (Kishimoto et al. 2011b; Burtscher et al. 2013).
On the other hand, the mid-IR interferometry for large angular size AGNs (i.e.
NGC1068 and Circinus galaxy) indicates an equatorial elongation in several
$R_{\rm sub}$ scales, embedded in the polar elongation at much larger scales
of tens of $R_{\rm sub}$ (Raban et al. 2009; Tristram et al. 2014; López-
Gonzaga et al. 2014). This has now even more clearly been shown by mid-IR
interferometric imaging for NGC1068 (Gámez Rosas et al. 2022; see more below)
and for Circinus (Isbell et al. 2022). Since we are probing a few $R_{\rm
sub}$ scale for NGC4151 here ($\lambda/b\sim 1.8$ mas; $R_{\rm sub}$
$\sim$0.3$-$0.5 mas for NGC4151), our detection of the equatorial elongation
seems geometrically consistent with such a structure seen by the mid-IR
interferometry.

In fact, with the near-IR interferometer VLTI/GRAVITY, Leftley et al. (2021)
has detected an equatorial elongation in the near-IR for the Type 1 AGN
ESO323-G77, despite the small visibility decrease at baselines up to 130 m.
The object shows a polar elongation in the mid-IR (Leftley et al., 2018). Thus
the equatorial elongation as seen by the CHARA Array for the type 1 AGN
NGC4151 seems to fit in with all these interferometric results.

The polar-elongated morphology found in the mid-IR at the pc scale, combined
with the larger-scale HST-resolved outflow, would suggest a conical and hollow
outflowing structure also at the inner pc scale. This wind has been suggested
to be driven by the radiation pressure on dust grains (Hönig et al., 2012,
2013). In this case, it is likely to be launched at the dust sublimation
region. There, the dusty gas could be inflated or puffed up by the infrared
radiation of the heated grains (Krolik, 2007), then radially accelerated by
the UV radiation from the nucleus and form the hollow conical structure
(Hönig, 2019). This near-IR emitting, wind-launching region might seem a
distinct, separate component in the spectral energy distribution (e.g. Mor et
al. 2009; Kishimoto et al. 2011b; Hönig et al. 2012), and in particular, this
could correspond to the bump-like feature at around 3 $\mu$m seen in the
spectra of many Type 1 AGNs (Hönig & Kishimoto 2017). Regarding this near-IR
emitting region as a slightly extended disk and the mid-IR emitting component
as a dusty wind, the model of Hönig & Kishimoto (2017) has successfully
explained the near-IR and mid-IR observations of the type 1 AGNs NGC3783 and
ESO323-G77 including interferometry (Leftley et al. 2021).

However, the observed emission properties of the 100-pc-scale outflow, the pc-
scale polar elongation, and the sub-pc-scale equatorially-elongated structure,
might actually suggest that all these three originate from the same nuclear
wind. They are all rather continuously distributed on the plane of the
surface-brightness and color temperature (Figure 7 in Miyauchi & Kishimoto
2020). We can also estimate the infrared emissivity (i.e. surface brightness
in units of Planck function) of these components on this plane, and infer the
optical thickness. While the 100 pc scale outflow is UV-optically thin, the
pc-scale structure (showing the polar elongation in the mid-IR) seems UV-
optically thick. This means that this inner hollow outflow would provide the
obscuration of the nucleus in equatorial directions in the UV/optical, as also
argued by Hönig (2019). In this case, the whole outflow would be the obscuring
‘torus’, and we would not need to invoke the disk component for obscuration
described above. The near-IR emitting component would not be a distinct one,
but would simply be the innermost part of the outflow. But then, we would
rather expect that it would look more polar-elongated, as opposed to the CHARA
Array results presented here.

###  5.5 Flaring geometry from episodic, one-off acceleration

One consistent picture might be that the distribution actually is
geometrically flaring effectively (see Figure 5a), with the cloud acceleration
being one-off, or episodic, rather than continuous. The geometrical flaring
would give an inner-equatorial, outer-polar structure, explaining the mid-IR
interferometry results for NGC1068 and Circinus described above with a single
structure (i.e. explains the two components of the inner-equatorial and outer-
polar structures as a single entity) and the CHARA Array results at the same
time. Here, we do not infer that each cloud follows the flaring path, as it
would need an additional acceleration toward the polar-axis direction in the
outer radii. Rather, in fact, it has been argued that the velocities of the
outflowing clouds at 100 pc scales could be from a one-off, episodic
acceleration (once in $\gtrsim$105 yr) at the innermost wind-launching site,
since they are Hubble-flow-like, i.e. the velocity of each cloud seems
proportional to its distance to the nucleus (Ozaki 2009; Miyauchi & Kishimoto
2020). The geometrical flaring above can realize if each cloud is simply
moving radially roughly at a constant speed but the clouds have a polar-fast,
equatorially-slow velocity distribution from the one-off anisotropic
acceleration (see Figure 5b), which might be reasonable for the radiation
pressure from an accretion disk.

This possible episodic acceleration might actually be quite consistent with
the episodic activity that has been inferred from X-ray observations. Wang et
al. (2010) has identified a significant soft X-ray emission extended out to
$\sim$$10^{4}$ light years from the nucleus of NGC4151. This emission has been
interpreted to indicate either mechanical heating by an episodic outflow or
photoionization by a highly-luminous, Eddington-limited outburst. In both
cases, the episodic activity is estimated to have occurred in the last
$10^{4}$-$10^{5}$ years, which is a timescale quite similar to what we discuss
here.

The flaring geometry from an episodic acceleration (Figure 5a) as seen edge-on
might fit with the mid-IR interferometric image in Gámez Rosas et al. (2022)
(their Figure 1e, mid-IR RGB-composite image). The illustration in Figure 5a,
which is from Figure 8 in Miyauchi & Kishimoto (2020), intends to represent
the distribution of the ionized gas and $T\gtrsim 300$ K dust, and the inner
equatorial region (represented in red) would correspond to the near-IR
emitting, dust sublimation region. This might also be consistent with the 3.7
$\mu$m image in Gámez Rosas et al. (2022) showing an equatorial structure,
with the dust sublimation region being still hidden below. In the case of the
Type 1 object NGC4151, we would see the sublimation region directly — Figure
5c shows the same 3D illustration, but rotated in 3D to match the inferred
inclination and PA of the axis projected on the sky for NGC4151 with the
central square schematically corresponding to the CHARA Array and KI
observations (Figure 5e). This illustrates the possible configuration in
NGC4151. Figure 5d shows the same geometry further rotated to have a view
perpendicular to both our line of sight and the system axis, corresponding to
the view of the reconstructed 3D distribution at a much larger scale shown in
Figure 4c and reproduced in Figure 5f.

###  5.6 Relations to other observations and simulations

How would this structure be related to those observed at other wavelengths and
other scales? As mentioned in section 5.1, the radio continuum source D3,
believed to contain the nucleus (Mundell et al. 2003), has been resolved at
$\sim$1 mas resolution with VLBI (a beam size of 0.95 $\times$ 0.46 mas;
Ulvestad et al. 2005). The source is resolved into three components along the
E-W radio jet direction. The middle component D3b is argued to be the nucleus,
with the other two knots $\sim$1$-$2 mas away from it. Thus the dust
sublimation region of $\sim$0.5 mas in radius, as resolved by the CHARA Array
observations here, would be located at the center of the three-knot jet-like
structure that has a scale comparable to that of the sublimation region
itself. This is schematically shown in Figure 5d, where approximate 3D
positions of the three knots are indicated. Here the jet axis at this scale is
assumed to coincide with the system axis, based on the inclination of
$\sim$40° inferred for the radio jet at larger scales (Pedlar et al. 1993)
which is consistent with the inclination of the system discussed in section
5.3.

At slightly larger, 10–100 mas scales, neutral hydrogen absorption is observed
in one of the several components of the linear radio jet. This is the
component called E+F, which is at the counter-jet side but closest to the
radio nucleus D3 (the projected distance is $\sim$70 mas or $\sim$5 pc;
Mundell et al. 1995, 2003). In Figure 5f, we have drawn approximate 3D
positions of the jet components assuming again the inclination of 40°. In
terms of the episodic polar-fast outflow discussed above, this foreground
neutral hydrogen material in front of the counter jet component E+F might be a
part of the slow-velocity outflows moving relatively close to equatorial
directions. In the scenario of the one-off acceleration $\sim$105 years ago,
the clouds several pc away from the nucleus would be at the velocity of
several 10s of km/s, which would be consistent with the observed radial
velocities of the absorption lines (Figure 5b of Mundell et al. 2003).

This central region seems to be surrounded by a molecular ring-like structure
at slightly outer radii of $\sim$1 arcsec, which is observed in the infrared
molecular hydrogen line H2 1-0 S(1), with the projected major axis PA
$\sim$170° and inclined by $\sim$37$\arcdeg$ from our line of sight (Fernandez
et al. 1999). This configuration is schematically shown in 3D in Figure 5f as
a red transparent disk (inner and outer radius of 0.5 and 1.5 arcsec,
respectively). The kinematics of this ring is inferred to be dominated by
rotation, rather than outflow, based on their Fabry-Perot image missing the
North-South part which would possess higher radial velocities if in rotation.
Therefore, this ring might be kinematically distinct from the obscuring and
outflowing torus discussed above, and rather be an outer part of the ‘fueling
disk’, an equatorial thin disk-like structure that transfers material toward
the nucleus (Hönig 2019; Miyauchi & Kishimoto 2020).

The idea to have a wind responsible for the obscuration of the nucleus has
been investigated with theoretical calculations and numerical simulations. The
obscuring wind could be from a combination of X-ray heating and radiation
pressure (Krolik & Begelman 1986, 1988; Wada 2012), or, as mentioned in
section 5.4, triggered by IR radiation pressure from heated dust and
accelerated by the AGN UV radiation (Krolik 2007; Chan & Krolik 2016, 2017;
Dorodnitsyn et al. 2016; Namekata & Umemura 2016; Williamson et al. 2019,
2020). As for the episodic property, the hydrodynamical winds simulated by
Wada (2012, 2015) are in fact quite intermittent in nature. Here we argued
observationally that an episodic wind accelerated by polar-strong radiation,
forming an effectively flaring geometry, would be able to explain the near-IR
equatorial elongation at a $R_{\rm sub}$ scale, mid-IR polar elongation at
tens of $R_{\rm sub}$ scale, as well as the linear velocity field and the
hollow conical outflow at a 100 pc scale, at least qualitatively. While the
exact physical processes for launching winds are still being explored
theoretically, a modeling focused on reproducing these observed properties
including possible episodic acceleration would enable us to do more
quantitative comparisons with observations. We plan to implement such modeling
as well as acquire further data on this innermost structure with the CHARA
Array.

##  6 Summary and conclusions

We have presented the near-IR interferometric measurements of the Type 1 AGN
NGC4151 at projected baselines up to $\sim$250 m, the longest baseline
achieved so far in the infrared for an extragalactic object. With these
baselines, the compact near-IR structure, or the dust sublimation region, has
now been spatially resolved. The measurements provide the size along different
position angles, giving the first information on the morphology of the object
at this inner spatial scale. When combined with the previous Keck
interferometer data albeit at shorter baselines, we find that the region is
elongated along roughly north-south, perpendicular to the system axis
direction as determined by the linear radio jet structure and optical
polarization observations.

We argued that an effectively flaring geometry can explain both the inner-
equatorial and outer-polar elongation simultaneously, which are seen in our
CHARA Array observations and the mid-IR interferometry. The flaring geometry
at a given time could be formed from a polar-fast, equatorially-slow, one-off
acceleration, while the episodic acceleration is quite consistent with the
Hubble-flow-like velocity field observed at a $\sim$100 pc scale. Long-
baseline measurements covering a wide position-angle range in a single epoch
are now desired to confirm the results. We are now aiming for acquiring such
data at the CHARA Array, which is currently the only infrared interferometer
that has such long baselines and sufficient uv coverage.

This work is based upon observations obtained with the Georgia State
University Center for High Angular Resolution Astronomy Array at Mount Wilson
Observatory. The CHARA Array is supported by the National Science Foundation
under Grant No. AST-1636624 and AST-2034336. Institutional support has been
provided from the GSU College of Arts and Sciences and the GSU Office of the
Vice President for Research and Economic Development. We are grateful for all
the generous and patient support from the CHARA staffs since the start of this
project in 2010 at the CHARA Array. A part of the time at the CHARA Array was
granted through the NOIRLab community access program (NOIRLab PropIDs:
2010A-0081, 2012A-0187, 2013A-0241, 2014A-0072, 2015A-0242, 2016A-0091,
2017A-0049, 2018B-0109, 2019A-0055, 2021A-0013; PI: M. Kishimoto). This work
was supported in part by JSPS grants 16H05731, 20K04029 and 21H04496, and the
grant E1906 through Kyoto Sangyo University. We acknowledge with thanks the
variable star observations from the AAVSO International Database contributed
by observers world-wide and used in this research. This research has made use
of the Jean-Marie Mariotti Center Aspro and SearchCal services.

## References

  * Antonucci (1983) Antonucci, R. R. J. 1983, Nature, 303, 158
  * Antonucci & Miller (1985) Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 297, 621, doi: 10.1086/163559
  * Anugu et al. (2020) Anugu, N., ten Brummelaar, T., Turner, N. H., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11446, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1144622
  * Burtscher et al. (2013) Burtscher, L., Meisenheimer, K., Tristram, K. R. W., et al. 2013, A&A, 558, A149, doi: 10.1051/0004-6361/201321890
  * Chan & Krolik (2016) Chan, C.-H., & Krolik, J. H. 2016, ApJ, 825, 67, doi: 10.3847/0004-637X/825/1/67
  * Chan & Krolik (2017) —. 2017, ApJ, 843, 58, doi: 10.3847/1538-4357/aa76e4
  * Che et al. (2013) Che, X., Sturmann, L., Monnier, J. D., et al. 2013, Journal of Astronomical Instrumentation, 2, 1340007, doi: 10.1142/S2251171713400072
  * Das et al. (2005) Das, V., Crenshaw, D. M., Hutchings, J. B., et al. 2005, AJ, 130, 945, doi: 10.1086/432255
  * Dorodnitsyn et al. (2016) Dorodnitsyn, A., Kallman, T., & Proga, D. 2016, ApJ, 819, 115, doi: 10.3847/0004-637X/819/2/115
  * Fernandez et al. (1999) Fernandez, B. R., Holloway, A. J., Meaburn, J., Pedlar, A., & Mundell, C. G. 1999, MNRAS, 305, 319, doi: 10.1046/j.1365-8711.1999.02406.x
  * Gámez Rosas et al. (2022) Gámez Rosas, V., Isbell, J. W., Jaffe, W., et al. 2022, Nature, 602, 403, doi: 10.1038/s41586-021-04311-7
  * GRAVITY Collaboration et al. (2019) GRAVITY Collaboration, Dexter, J., Shangguan, J., et al. 2019, arXiv e-prints, arXiv:1910.00593.  https://arxiv.org/abs/1910.00593
  * Hönig (2019) Hönig, S. F. 2019, ApJ, 884, 171, doi: 10.3847/1538-4357/ab4591
  * Hönig & Kishimoto (2017) Hönig, S. F., & Kishimoto, M. 2017, ApJ, 838, L20, doi: 10.3847/2041-8213/aa6838
  * Hönig et al. (2012) Hönig, S. F., Kishimoto, M., Antonucci, R., et al. 2012, ApJ, 755, 149, doi: 10.1088/0004-637X/755/2/149
  * Hönig et al. (2013) Hönig, S. F., Kishimoto, M., Tristram, K. R. W., et al. 2013, ApJ, 771, 87, doi: 10.1088/0004-637X/771/2/87
  * Isbell et al. (2022) Isbell, J. W., Meisenheimer, K., Pott, J. U., et al. 2022, A&A, 663, A35, doi: 10.1051/0004-6361/202243271
  * Kishimoto et al. (2011a) Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011a, A&A, 527, A121, doi: 10.1051/0004-6361/201016054
  * Kishimoto et al. (2009) —. 2009, A&A, 507, L57.  https://arxiv.org/abs/0911.0666
  * Kishimoto et al. (2011b) —. 2011b, A&A, 536, A78, doi: 10.1051/0004-6361/201117367
  * Kishimoto et al. (2007) Kishimoto, M., Hönig, S. F., Beckert, T., & Weigelt, G. 2007, A&A, 476, 713, doi: 10.1051/0004-6361:20077911
  * Kishimoto et al. (2013) Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2013, ApJ, 775, L36, doi: 10.1088/2041-8205/775/2/L36
  * Koshida et al. (2009) Koshida, S., Yoshii, Y., Kobayashi, Y., et al. 2009, ApJ, 700, L109, doi: 10.1088/0004-637X/700/2/L109
  * Krolik (2007) Krolik, J. H. 2007, ApJ, 661, 52, doi: 10.1086/515432
  * Krolik & Begelman (1986) Krolik, J. H., & Begelman, M. C. 1986, ApJ, 308, L55, doi: 10.1086/184743
  * Krolik & Begelman (1988) —. 1988, ApJ, 329, 702, doi: 10.1086/166414
  * Leftley et al. (2021) Leftley, J. H., Tristram, K. R. W., Hönig, S. F., et al. 2021, ApJ, 912, 96, doi: 10.3847/1538-4357/abee80
  * Leftley et al. (2018) —. 2018, ApJ, 862, 17, doi: 10.3847/1538-4357/aac8e5
  * Lira et al. (2020) Lira, P., Goosmann, R. W., Kishimoto, M., & Cartier, R. 2020, MNRAS, 491, 1, doi: 10.1093/mnras/stz2774
  * López-Gonzaga et al. (2016) López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, A&A, 591, A47, doi: 10.1051/0004-6361/201527590
  * López-Gonzaga et al. (2014) López-Gonzaga, N., Jaffe, W., Burtscher, L., Tristram, K. R. W., & Meisenheimer, K. 2014, A&A, 565, A71, doi: 10.1051/0004-6361/201323002
  * Martel (1998) Martel, A. R. 1998, ApJ, 508, 657, doi: 10.1086/306453
  * Minezaki et al. (2004) Minezaki, T., Yoshii, Y., Kobayashi, Y., et al. 2004, ApJ, 600, L35, doi: 10.1086/381364
  * Miyauchi & Kishimoto (2020) Miyauchi, R., & Kishimoto, M. 2020, ApJ, 904, 149, doi: 10.3847/1538-4357/abc264
  * Mor et al. (2009) Mor, R., Netzer, H., & Elitzur, M. 2009, ApJ, 705, 298, doi: 10.1088/0004-637X/705/1/298
  * Mundell et al. (1995) Mundell, C. G., Pedlar, A., Baum, S. A., et al. 1995, MNRAS, 272, 355, doi: 10.1093/mnras/272.2.355
  * Mundell et al. (2003) Mundell, C. G., Wrobel, J. M., Pedlar, A., & Gallimore, J. F. 2003, ApJ, 583, 192, doi: 10.1086/345356
  * Namekata & Umemura (2016) Namekata, D., & Umemura, M. 2016, MNRAS, 460, 980, doi: 10.1093/mnras/stw862
  * Ozaki (2009) Ozaki, S. 2009, PASJ, 61, 259, doi: 10.1093/pasj/61.2.259
  * Pedlar et al. (1993) Pedlar, A., Kukula, M. J., Longley, D. P. T., et al. 1993, MNRAS, 263, 471, doi: 10.1093/mnras/263.2.471
  * Pott et al. (2010) Pott, J., Malkan, M. A., Elitzur, M., et al. 2010, ApJ, 715, 736, doi: 10.1088/0004-637X/715/2/736
  * Raban et al. (2009) Raban, D., Jaffe, W., Röttgering, H., Meisenheimer, K., & Tristram, K. R. W. 2009, MNRAS, 394, 1325, doi: 10.1111/j.1365-2966.2009.14439.x
  * Schnülle et al. (2015) Schnülle, K., Pott, J. U., Rix, H. W., et al. 2015, A&A, 578, A57, doi: 10.1051/0004-6361/201525733
  * Smith et al. (2004) Smith, J. E., Robinson, A., Alexander, D. M., et al. 2004, MNRAS, 350, 140, doi: 10.1111/j.1365-2966.2004.07610.x
  * Swain et al. (2003) Swain, M., Vasisht, G., Akeson, R., et al. 2003, ApJ, 596, L163, doi: 10.1086/379235
  * ten Brummelaar et al. (2005) ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453, doi: 10.1086/430729
  * Ten Brummelaar et al. (2013) Ten Brummelaar, T. A., Sturmann, J., Ridgway, S. T., et al. 2013, Journal of Astronomical Instrumentation, 2, 1340004, doi: 10.1142/S2251171713400047
  * ten Brummelaar et al. (2018) ten Brummelaar, T. A., Sturmann, J., Sturmann, L., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10703, Adaptive Optics Systems VI, ed. L. M. Close, L. Schreiber, & D. Schmidt, 1070304
  * Tristram et al. (2014) Tristram, K. R. W., Burtscher, L., Jaffe, W., et al. 2014, A&A, 563, A82, doi: 10.1051/0004-6361/201322698
  * Ulvestad et al. (2005) Ulvestad, J. S., Wong, D. S., Taylor, G. B., Gallimore, J. F., & Mundell, C. G. 2005, AJ, 130, 936, doi: 10.1086/432034
  * Wada (2012) Wada, K. 2012, ApJ, 758, 66, doi: 10.1088/0004-637X/758/1/66
  * Wada (2015) —. 2015, ApJ, 812, 82, doi: 10.1088/0004-637X/812/1/82
  * Wang et al. (2010) Wang, J., Fabbiano, G., Risaliti, G., et al. 2010, ApJ, 719, L208, doi: 10.1088/2041-8205/719/2/L208
  * Weigelt et al. (2012) Weigelt, G., Hofmann, K.-H., Kishimoto, M., et al. 2012, A&A, 541, L9, doi: 10.1051/0004-6361/201219213
  * Williamson et al. (2019) Williamson, D., Hönig, S., & Venanzi, M. 2019, ApJ, 876, 137, doi: 10.3847/1538-4357/ab17d5
  * Williamson et al. (2020) —. 2020, ApJ, 897, 26, doi: 10.3847/1538-4357/ab989e

##  Appendix A Observed raw visibilities and power spectra of fringe tracks

Figure 6: Observed raw visibilities in the night of 2020-02-15 at CHARA. Light
dotted lines show the raw visibilities averaged over 10 scans ($\sim$15
seconds) in each track of objects. Blue filled circles show the raw
visibilities integrated over each track of calibrators with uncertainties.
Blue open circles are the system visibilities estimated from the calibrator
observations. The green filled circle is the raw visibility of the science
target from the integration over the track.

Figure 7: Details of the fringe track for NGC4151 on 2020-02-15. (a) Fringe
envelopes of the scans are shown in grayscale, with each scan in horizontal
direction, forming a waterfall plot along the vertical direction. (b) Power
spectrum of each scan is shown in grayscale along the horizontal direction,
again forming a waterfall plot of power spectra. (c) Dotted lines show the raw
$V^{2}$ obtained from integrating each power spectrum, and the blue dots
indicate its 10-scan average. (d) Averaged power spectrum is shown in red
dotted line, with its smoothed version in red solid line. The blue spectrum is
before the subtraction of the noise bias, which is shown in gray, with the
scale axis at the right side in the same blue.

Figure 7 shows the raw visibilities for each calibrator and the target
observed in the night of 2020-02-015. For each object, the raw visibilities
obtained by averaging the power spectra over 10 scans (corresponding to an
integration over $\sim$15 seconds) are shown in dotted lines. The visibilities
integrated over each track of calibrator observations are shown in filled blue
circles with uncertainties. The open circles show the system visibilities
(transfer functions) at the time of the calibrator observations (almost
overlapping with the calibrators’ raw visibilities) and at the time of the
science target observation. The latter is obtained from the linear
interpolation between the system visibilities bracketing the science target
observation. The visibility integrated over the track (or a portion of the
track; see below) for the science target is shown in a green filled circle
with its uncertainty.

Figure 7 shows the details of the fringe track for NGC4151 on 2020-02-15. The
top-left panel (a) shows the envelope (fringe amplitude) over each scan in
grayscale (note the right-side axis showing approximate elapsed time in
seconds). The top-middle panel (b) shows the corresponding power spectra also
in grayscale. Fringes should be detected at the frequency of 1 cycle per
observing wavelength $\lambda$. They are clearly seen with some wiggles around
this frequency over different scans. The top-right panel (c) shows the raw
$V^{2}$ value for each scan in dotted lines, while the light blue dots
correspond to the averages over 10 scans. In the bottom panel (d), the red
dotted line is the final averaged power spectrum, with the same spectrum
smoothed over 10 frequency pixels shown in red solid line. This is after the
subtraction of the noise bias spectrum, which is measured from the average of
the same number of off-fringe scans as the object scans. Here the power
spectrum before the subtraction is shown in blue and the noise bias spectrum
in gray, with the scale axis for both indicated at the right side in blue.

In this particular fringe track, the fringes slipped away from the scan range
at and around the scan #220. Therefore, we integrated the power spectra until
scan #190. The integration range is indicated with a green vertical line in
panel (c), with its x-axis position corresponding to the value of the obtained
$V^{2}$ from the integration. The uncertainty is also shown in the horizontal
green line at the mid-point of integration. Note that the exact choice of the
scan number for integration does not alter the results beyond the estimated
uncertainty.

Figures 9-16 show the observed raw visibilities in the same way. Multiple
tracks of the science target within a given night are marked with sequential
numbers, and shown individually in separate Figures. Note that we did not see
fringes in the track #01 for NGC4151 in the night of 2021-04-30. Since we are
not entirely sure if the scans were in the right range, this track was
excluded from the analysis, but shown here for completeness.

Figure 8: Observed raw visibilities in the night of 2020-02-16 at CHARA. The
notations are the same as in Figure 7.

Figure 9: Details of the fringe track for NGC4151 on 2020-02-16, in the same
format as in Figure 7.

Figure 10: Observed raw visibilities in the night of 2021-03-19 at CHARA. The
notations are the same as in Figure 7.

Figure 11: Details of the fringe track for NGC4151 on 2021-03-19, in the same
format as in Figure 7. Figure 12: Observed raw visibilities in the night of
2021-04-30 at CHARA. The notations are the same as in Figure 7. Figure 13:
Details of the fringe track #01 for NGC4151 on 2021-04-30, in the same format
as in Figure 7. Figure 14: Details of the fringe track #02 for NGC4151 on
2021-04-30, in the same format as in Figure 7. Figure 15: Details of the
fringe track #03 for NGC4151 on 2021-04-30, in the same format as in Figure 7.
Figure 16: Details of the fringe track #04 for NGC4151 on 2021-04-30, in the
same format as in Figure 7.