# TaylorImNet for Fast 3D Shape Reconstruction Based on
Implicit Surface Function

Yuting Xiao1 Jiale Xu111footnotemark: 1 Shenghua Gao1,2
1ShanghaiTech University
2Shanghai Engineering Research Center of Intelligent Vision and Imaging
{xiaoyt, xujl1<EMAIL_ADDRESS>
Equal contribution.Corresponding author.

###### Abstract

Benefiting from the contiguous representation ability, deep implicit functions
can extract the iso-surface of a shape at arbitrary resolution. However,
utilizing the neural network with a large number of parameters as the implicit
function prevents the generation speed of high-resolution topology because it
needs to forward a large number of query points into the network. In this
work, we propose TaylorImNet inspired by the Taylor series for implicit 3D
shape representation. TaylorImNet exploits a set of discrete expansion points
and corresponding Taylor series to model a contiguous implicit shape field.
After the expansion points and corresponding coefficients are obtained, our
model only needs to calculate the Taylor series to evaluate each point and the
number of expansion points is independent of the generating resolution. Based
on this representation, our TaylorImNet can achieve a significantly faster
generation speed than other baselines. We evaluate our approach on
reconstruction tasks from various types of input, and the experimental results
demonstrate that our approach can get slightly better performance than
existing state-of-the-art baselines while improving the inference speed with a
large margin. Our code will be released in https://github.com/svip-
lab/TaylorImNet.

##  1 Introduction

Figure 1: (a) The illustration of our Taylor implicit functions. We represent
a shape with a set of expansion points (Hollow Black) and corresponding Taylor
series. Each expansion point controls a local region around it. The signed
distance of a query point $\boldsymbol{x}$ is computed by the weighted average
of the Taylor series at the nearby expansion points $A,B$ and $C$. (b) The
evaluation of the representation capacity of the Taylor series. We separately
sample 1000 expansion points from 1300 objects in the ShapeNet dataset and
apply the least-squares method to estimate the Taylor series coefficients of
each expansion point. We visualize the mean error of SDF estimated by
different-order Taylor series, the results demonstrate that low-order Taylor
series can already represent the shape accurately.

Deep implicit functions parameterized by neural networks have gained much
popularity in 3D shape reconstruction and representation. Compared with
explicit shape representations such as voxels [14, 17, 38], point clouds [30,
31, 1, 41] and meshes [37, 22], deep implicit functions represent the shape
surface as the continuous decision boundary. At the inference phase, deep
implicit functions predict the occupancy or signed distance field (SDF) for
each querying position and extract the iso-surface with a post-processing
step, _i.e_., marching cubes [24]. The implicit functions separate each point
in the space into inside and outside and can generate a continuous iso-surface
for each shape. Some previous works [25, 6, 33] apply neural network as the
implicit function to generalize all shapes, while some other works [19, 13]
observe that the global 3D implicit representation may be hard to model and
generalize because the possible global shapes are too diverse and might be too
complex. So these works attempt to use the deep network to infer the local 3D
representation because the local appearances of different objects are similar
and easy to generalize.

However, deep implicit functions need to forward the neural network for every
query point at the inference phase, which is computation-unfriendly when
generating surface at high resolution since there will be a large amount of
points to query. Enhancing the capacity of the neural network to model complex
shapes inevitably leads to increasing computational requirements and it is
difficult to balance the capacity demand of the neural network and inference
speed. Besides, though some works have discussed using a deep network with a
large capacity to infer the parameters of a simple network for modeling the
local region, which is similar to the conditional convolution, it is still
tough to apply these approaches in reconstruction because the redundancy of
parameters makes the parameter prediction too difficult, even it is a very
simple deep network.

In this paper, inspired by the Taylor expansion, we propose a Taylor-series-
based compact implicit local 3D shape representation to solve the
aforementioned problem. Each object is represented by an expansion point set
and the corresponding Taylor series. Furthermore, we propose a coarse-to-fine
expansion points sampling strategy in the inference phase. The number of
expansion points is independent of the resolution of mesh generation. Instead
of utilizing the neural network to evaluate every query point, adopting the
discrete point set and corresponding Taylor series to represent 3D shapes
makes the inference independent to the neural network once the expansion
points and corresponding coefficients are obtained. This representation can
also take advantage of the learning ability of the neural network meanwhile
promoting the computational efficiency when generating high-resolution
topology since the time of inferring the coefficients of the Taylor series is
fixed for each object and it does not need to forward the network for new
query points.

Applying a single expansion point to cover a local region makes the global
shape field discrete. To model the contiguous surface, we compute the weighted
average of the Taylor series among several nearest expansion points of each
query point, as shown in Figure 1 (a). The union of local regions around every
expansion point can cover the global shape region. After the expansion points
and the corresponding Taylor series coefficients are obtained, we can infer
the signed distance field of arbitrary points in the space.

We evaluate our approach on the 3D object reconstruction task with 3 different
input types: single-view images, point clouds, and voxels. In our experiments,
we observe that only applying low-order Taylor series can reconstruct shape
accurately, which is consistent with the experimental results in Figure 1 (b).
And the inference speed is improved with a large margin compared with
classical implicit function-based baselines.

The contributions of this paper could be summarized as:

  * •

We propose a compact and computational-friendly implicit 3D shape
representation based on the Taylor expansion.

  * •

We propose a coarse-to-fine expansion points sampling strategy in the
inference phase. This process is independent of the resolution of mesh
generation.

  * •

We demonstrate our implicit 3D representation can significantly speed up the
generation while generating meshes with slightly better quality than existing
state-of-the-art baselines on different input data types.

##  2 Related Work

3D Shape Reconstruction. There has been significant work on 3D shape
reconstruction from various types of shape observations. Given a single-view
image, 3D-R2N2 [8], PSGN [11], AtlasNet [16], Pixel2Mesh [37] and Pix2Vox [39]
can reconstruct a explicitly represented model, _e.g_., voxel grids, point
clouds and meshes. Compared to explicit shape representations, implicit
functions parameterized by neural networks is more memory-efficient and can
preserve fine-grained shape details, thus they are adopted by many methods,
_e.g_., OccNet [25], and PIFu [33]. Different from single-view reconstruction
which is an ill-posed problem, reconstruction from discrete point clouds or
coarse voxel grids leads to finer shapes due to the 3D prior knowledge
provided by them inherently. Given a point cloud as input, traditional
optimization-based methods such as Moving Least Square [2] and Poisson Surface
Reconstruction [20, 21], and deep optimization-based methods such as SAL [3]
and IGR [15], can generate a continuous surface with fine details. When ground
truth point-occupancy pairs are available, implicit-function-based methods
such as OccNet [25] and ConvOccNet [29] can also address this task well.
Besides, these implicit-function-based models [25, 29, 7] can also reconstruct
high-resolution surfaces from coarse voxel grids (_i.e_., voxel super-
resolution) by assembling a proper 3D convolutional encoder. Despite the
preferable performance, implicit-function-based methods need a post-processing
step (_e.g_., Marching Cubes [24]) to extract the iso-surface, which is often
time-consuming at high resolution.

Figure 2: The Overview of our model. The input data types of this framework
can be the single-view, the point cloud, or the voxel. In the training phase,
Each expansion point $\boldsymbol{x}^{\prime}_{i}$ has a set of sampled points
nearby and the signed distance of each point $\boldsymbol{x}_{i,j}$ in this
set is computed based on the implicit Taylor series at expansion point
$\boldsymbol{x}^{\prime}_{i}$. In the inference phase, after the expansion
points sampling, the signed distance of arbitrary point $\boldsymbol{x}$ are
inferred by the weighted average of the implicit Talor expansion functions
corresponding to its multiple nearest expansion points.

Global Shape Representations. Global shape representations seek to represent a
3D shape with a single implicit latent code. Following this spirit, an auto-
encoder (AE) architecture is often exploited [25, 6]. Specifically, an encoder
maps the input shape observation into the latent space, and a decoder tries to
recover the shape with the holistic latent code. In comparison, DeepSDF [27]
adopts an auto-decoder scheme, which randomly samples a latent code from a
Gaussian prior distribution for each shape and optimizes it with gradient
descent during the training process. A series of work [10, 26, 42, 9, 34] has
made improvements upon the framework of DeepSDF and has shown promising shape
modeling results. Global shape representations have good shape completion
ability and the learned latent space is convenient for shape interpolation and
generation. However, they struggle to recover fine details. Thus, some recent
work [35, 36] seek to use periodic activation functions for better preserving
high-frequency details.

Local Shape Representations. Local shape representations decompose a holistic
shape into local parts to better preserve fine details. Some methods learn to
decompose the shape into local parts automatically, and represent the parts
with CSG primitives [32], superquadrics [28] or 3D Gaussians [12, 13]. Dualsdf
[18] expresses shapes at two levels of granularity, one capturing fine details
and the other representing an abstracted proxy shape. Another set of work
seeks to divide the 3D space into local patches. DeepLS [4] uniformly divides
the 3D space into voxel grids and assigns a latent code as well as a DeepSDF
decoder to each voxel. LGCL [40] adopts a similar idea except that the space
partition is defined by a set of key points. LIG [19] train a part auto-
encoder to learn an embedding of local crops of 3D shapes, and then the
decoder is utilized to optimize the latent code of local crops at inference
time.

##  3 Method

In this section, we introduce the modeling of our Taylor-based implicit 3D
shape representation, how to train our approach, and how to utilize our
coarse-to-fine to sample expansion points in the inference.

The classical deep implicit function based reconstruction methods first learn
a latent representation $\boldsymbol{c}$ from the input partial 3D
observations (_e.g_., images, point clouds, voxels, and SDF grids), and then
utilize a neural network to predict $O$ (the SDF or occupancy value) of a
query point $\boldsymbol{x}$ from the latent representation. This can be
denoted as:

$O=f(\boldsymbol{x},\mathbf{c})$ (1)

Since that the $f$ need to forward once for every query point, it is
computationally unfriendly to generate high-resolution meshes.

Our approach aims to utilize a fixed set of expansion points and corresponding
Taylor series to represent a 3D shape. After these are obtained, the inferring
of points is independent of the network-based decoder and this significantly
speeds up the generation speed, especially for the high-resolution mesh.

###  3.1 Taylor-based Implicit Function Formulation

Given an signed distance field $F(\boldsymbol{x})$ of a shape, where the
$\boldsymbol{x}\in\mathbb{R}^{3}$ is the point coordinate. Our target is to
utilize the combination of the Taylor series to represent the approximation of
the $F(\boldsymbol{x})$. Suppose that
$\boldsymbol{x}^{\prime}_{i}\in\mathbb{R}^{3}$ is an expansion point. We apply
the Taylor series at this expansion point $\boldsymbol{x}^{\prime}_{i}$ to
model the target shape $F(\boldsymbol{x})$ at the local region near
$\boldsymbol{x}^{\prime}_{i}$. The 3D representation of local region near the
expansion point $\boldsymbol{x}^{\prime}_{i}$ is denoted as
$F_{s}(\boldsymbol{x},\boldsymbol{x}^{\prime}_{i})$:

$\displaystyle
T(\boldsymbol{x}^{\prime}_{i})=(F(\boldsymbol{x}^{\prime}_{i}),\nabla_{\boldsymbol{x}}F(\boldsymbol{x}^{\prime}_{i}),\frac{1}{2}\nabla_{\boldsymbol{x}}^{2}F(\boldsymbol{x}^{\prime}_{i}),...)$
(2) $\displaystyle
X(\boldsymbol{x},\boldsymbol{x}^{\prime}_{i})=(1,\boldsymbol{x}-\boldsymbol{x}^{\prime}_{i},(\boldsymbol{x}-\boldsymbol{x}^{\prime}_{i})(\boldsymbol{x}-\boldsymbol{x}^{\prime}_{i})^{T},...)$
(3) $\displaystyle
F_{s}(\boldsymbol{x};\boldsymbol{x}^{\prime}_{i})=T(\boldsymbol{x}^{\prime}_{i})\cdot
X(\boldsymbol{x},\boldsymbol{x}^{\prime}_{i})$ (4)

where the $T(\boldsymbol{x}^{\prime}_{i})$ is the coefficients of the Taylor
series at expansion point $\boldsymbol{x}^{\prime}_{i}$.

For an arbitrary point $\boldsymbol{x}$ in the space, inferring the
$F(\boldsymbol{x})$ from the single nearest expansion point would make the
signed distance field discontinuous. To alleviate the discontinuities, the
$F(\boldsymbol{x})$ can be defined by the weighted average of the signed
distance functions corresponding to multiple expansion points. Suppose that
there is a set of expansion points $\mathcal{P}$ close to the point
$\boldsymbol{x}$ in the space, the $F(\boldsymbol{x})$ can be noted as:

$F(\boldsymbol{x})={\sum}_{\boldsymbol{x}^{\prime}_{i}\in\mathcal{P}}\omega(\|\boldsymbol{x}-\boldsymbol{x}^{\prime}_{i}\|,\theta)F_{s}(\boldsymbol{x};\boldsymbol{x}^{\prime}_{i})$
(5)

where the $\mathcal{P}$ denotes the set of expansion points in the $n$ nearest
neighbourhood of point $\boldsymbol{x}$, and the weight function is set as the
Softmin function: $\omega(d_{i},\theta)=\text{exp}(-\theta
d_{i})/\sum_{i}\text{exp}(-\theta d_{i})$ , since the expansion point far from
the query point $x$ should obtain large weight decay. $d_{i}$ is the distance
from the point $\boldsymbol{x}$ to the nearby expansion point
$\boldsymbol{x}^{\prime}_{i}$.

###  3.2 Training Phase

We utilize an encoder $f_{e}$ and a decoder $f_{d}$ to infer the parameters of
our Taylor series based 3D shape representation.

We randomly sample $N$ expansion points in the 3D volume of an object. For
each expansion point $\boldsymbol{x}^{\prime}_{i}$ sampled in the space, we
input the corresponding feature code $\mathbf{E}$ and the position of this
expansion point to the decoder $f_{d}$, and the output is the parameters
$T(\boldsymbol{x}^{\prime}_{i})$ of our implicit shape representation, which
is denoted as:

$T(\boldsymbol{x}^{\prime}_{i})=f_{d}(\boldsymbol{c},\boldsymbol{x}^{\prime}_{i})$
(6)

After the parameters of an expansion point $\boldsymbol{x}^{\prime}_{i}$ are
obtained, we can predict the signed distance of any point $\boldsymbol{x}$
nearby this expansion point:

$s(\boldsymbol{x},\boldsymbol{x}^{\prime}_{i})=T(\boldsymbol{x}^{\prime}_{i})\cdot
X(\boldsymbol{x},\boldsymbol{x}^{\prime}_{i})^{T}$ (7)

Reconstruction loss. Directly applying the signed distance field for
supervising the learning of the network is not applicable because the signed
distance error at the point far from the surface is not significant for
reconstruction. Aiming to lead the network to focus on the reconstruction of
the surface region, we apply a sigmoid transformation on the signed distance
field:

$\sigma(s)=\frac{1}{1+e^{-\alpha s}}$ (8)

where $\alpha$ is a hyper-parameter. After the transformation, the points far
from the surface would have almost the same distance value, and the decline in
diversity makes the prediction easier than before the transformation. Besides,
applying the TSDF (Truncated Signed Distance Field) and binary occupancy field
has the same effect but both of them suffer from the existence of the non-
differentiable points, which needs high order of Taylor series to approximate.

To supervise the inferred Taylor series coefficients
$T(\boldsymbol{x}^{\prime}_{i})$, we sample a grid of points centered around
the expansion point $\boldsymbol{x}^{\prime}_{i}$ and generate the ground
truth signed distance values of these points. Suppose the
$\boldsymbol{x}_{i,j}\in\mathbb{R}^{3}$ is the $j$-th point in the point set
sampled in the cube of the $i$-th expansion point
$\boldsymbol{x}^{\prime}_{i}$ and $M$ points are sampled in the cube for each
expansion point. We design the reconstruction loss function in the training
phase as:

$\mathcal{L}=\frac{1}{NM}\sum_{i=0}^{N}\sum_{j=0}^{M}\mathcal{L}_{ce}(\sigma(s(\boldsymbol{x}_{i,j},\boldsymbol{x}^{\prime}_{i})),\sigma(F(\boldsymbol{x}_{i,j})))$
(9)

where the $N,M$ is the number of the expansion points and the sampled points
in each group, and the $\mathcal{L}_{ce}$ is the cross-entropy loss.

Figure 3: An example to introduce the coarse-to-fine strategy for expansion
point sampling in the inference phase. A point is estimated by 4 nearest
expansion points with the weights are $0.79$, $0.11$, $0.06$, and $0.04$
respectively.

###  3.3 Inference Phase

In this section, we introduce the coarse-to-fine strategy to obtain the iso-
surface of each object in the inference phase.

Firstly, a grid is generated to localize the coarse area of the surface. We
divide the 3D volume with a low-resolution grid for sampling expansion points.
For each unit of the grid, we sample the center point as the expansion point
$\boldsymbol{x}^{\prime}_{i}$ and compute the signed distance of this
expansion point. An inside distance threshold $\epsilon_{in}$ and an outside
threshold $\epsilon_{out}$ is set to filter the units into three classes: 1.
outside units far from the surface; 2. inside unit far from the surface; 3.
units near the surface. we aim to only sample dense expansion points in the
units close to the surface, which can be denoted as:

$\mathcal{R}=\\{R(\boldsymbol{x}^{\prime}_{i})|\epsilon_{in}\leq\sigma(s(\boldsymbol{x}^{\prime}_{i},\boldsymbol{x}^{\prime}_{i}))\leq\epsilon_{out}\\}$
(10)

where the $R(\boldsymbol{x}^{\prime}_{i})$ is the region of the unit
corresponding to the expansion point $\boldsymbol{x}^{\prime}_{i}$. The
$s(\boldsymbol{x}^{\prime}_{i},\boldsymbol{x}^{\prime}_{i})$ can be consider
as the Taylor expansion with 0 order at $\boldsymbol{x}^{\prime}_{i}$

Then, after the coarse surface regions are obtained, we sample more expansion
points in each unit near the surface. This operation utilizes more expansion
points to make the surface reconstruction more fine and smooth.

As shown in Figure 3. Based on the implicit 3D representation inferred by our
model, we can infer the signed distance of arbitrary points in the space. For
the point belonging to a unit far from the surface, the occupancy of this
point is set to the same as the unit it belongs to. For the point belonging to
the unit near the surface, we choose its $k$ nearest neighborhood expansion
points to compute the weighted average signed distance. The Marching Cubes
algorithm to obtain the predicted iso-surface.

##  4 Experiments

In this section, we conduct experiments to validate the performance of our
TaylorImNet. Firstly, we compare the reconstruction results on 3 different
types of input quantitatively and qualitatively. Secondly, we demonstrate the
inference acceleration effect compared with classical implicit-function-based
methods. Thirdly, we conduct ablation studies to evaluate the representation
capacity of our approach.

Dataset. Our experiments are conducted on the ShapeNet [5] dataset which
contains 13 shape categories. We follow the data pre-processing operations in
OccNet [25] for the fair comparison. For the single-view reconstruction and
the voxel super-resolution task, we use the image renderings and model
voxelizations provided by 3D-R2N2 [8]. For the point cloud reconstruction
task, the point clouds sampled by [25] are exploited as input.

Metrics. Following the ConvOccNet[29], we quantify the performance of our
method and other baselines with three metrics: (i) volumetric IoU (higher is
better), (ii) Chamfer-$L_{1}$ distance (lower is better), and (iii) F-Score
(higher is better). To compute the volumetric IoU, 100k points are uniformly
sampled in the bounding volume of the mesh. The Chamfer-$L_{1}$ distance is
computed by randomly sampling 100k points on the mesh surface and we scale the
Chamfer-$L_{1}$ results by 10. The F-Score is reported with threshold
$\tau=0.01$.

Training Details. In the training phase, we sample 2 types of points: 1.
expansion points for shape representation; 2. query points near each expansion
point for supervising the prediction of Taylor series coefficients in Equation
9. For the expansion point sampling, we uniformly random sample $k_{1}$
expansion points in the bounding volume and $k_{2}$ expansion points near the
surface. We set $k_{1}:k_{2}=1:3$ in our implementation. For the query point
sampling, we aim to sample a set of points near each expansion point. Treating
a chosen expansion point as the center point, we evenly sample a $5^{3}$ grid
of query points in a cube whose side length equals 0.08. Based on this, each
expansion point has 125 query points nearby. Besides, the $\alpha$ of the
sigmoid transformation is set to 32. We use the Adam optimizer with
$\beta_{1}=0.9$ and $\beta_{2}=0.999$. We train our model with a learning rate
starting from $1\times 10^{-3}$ and divided by 10 for 2 times in the training.
The implementation details of different data types are introduced in the
supplementary material.

Inference Details. In the inference phase, take the single-view reconstruction
as an example, our coarse-to-fine expansion point sampling strategy firstly
samples a $16^{3}$ grid of expansion points and we set $\epsilon_{out}=0.98$
and $\epsilon_{in}=0.02$ to filter the units near the surface. Then, we sample
a $2^{3}$ grid in each unit near the surface. After the expansion points are
obtained, we apply the 4 nearest expansion points to evaluate the point in the
unit near the surface.

|  IoU $\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score $\uparrow$
---|---|---|---
3D-R2N2[8] | 0.500 | 0.246 | 0.347
PSGN[11] | - | 0.215 | 0.142
Pixel2Mesh [37] | 0.480 | 0.216 | 0.274
OccNet [25] | 0.593 | 0.194 | 0.523
Ours | 0.599 | 0.192 | 0.553
Table 1: The quantitative results of the single-view reconstruction task.
Figure 4: Qualitative results on single-view reconstruction task.

###  4.1 Reconstruction from Single-view Images

We first conduct the single-view reconstruction experiments. The image
renderings provided by [8] are used as input, which contains 24 views for each
shape. We compare our model with 4 baseline approaches: 3D-R2N2 [8], PSGN
[11], Pixel2Mesh [37] and OccNet [25]. Our model uses the same backbone as
OccNet and achieve slightly better performance, which indicates the efficient
capacity of shape representation. The quantitative results are shown in Table
1. Our approach can obtain comparable performance on IoU and Chamfer-$L_{1}$
metrics.

We also show some qualitative results in Figure 4. We observe that our model
achieves better performance on curved surface reconstruction than other
compared baselines, such as the curved table legs and sleek chair armrests as
shown in the 2nd and 3rd rows respectively. This benefited from the high-order
terms of the Taylor series which can easily represent the curved surface
compared with OccNet [25]. The pixel2mesh [37] can create sharp corners while
hardly reconstructing smooth curved surfaces. And the table legs in the first
row show that our approach can reconstruct reasonable results on straight and
sharp parts. The 4th and 5th rows show that our approach has enough ability to
reconstruct the right angle area such as the right-angle handrail and pillar
of chairs.

###  4.2 Reconstruction from Point Clouds

This task aims to reconstruct the mesh surface from a point cloud with noise.
We randomly subsample 3000 points from point cloud provided by OccNet [25] as
input and apply Gaussian noise with zero mean and standard deviation 0.005.
The implementation of DMC [23] and PSGN [11] is adopted from the OccNet
codebase111https://github.com/autonomousvision/occupancy_networks. Our model
uses the same backbone as ConvOccNet [29] with its volume encoder.

Figure 5: Qualitative results on point cloud reconstruction task. The red
dotted boxes highlight the advantages of our approach.

The quantitative results are shown in Figure 5. Our approach and ConvOccNet
can obtain obviously more reasonable results than other baselines. And
compared with ConvOccNet, our approach has a significant advantage for
generating flat and narrow shape parts and holes. The red dotted boxes in the
last 3 columns demonstrate that ConvOccNet can not reconstruct reasonable
surfaces at the local region where the occupancy of points changes frequently.
This situation happens when the different parts of a surface are stacked with
each other. For example, the backplane of the tail of the airplane (1st row),
the mezzanine under a chair cushion or the table(2nd, 4th, and 5th rows), and
the region between the back of the display and the pillar (3rd row).

|  IoU $\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score $\uparrow$
---|---|---|---
PSGN [11] | - | 0.178 | 0.180
DMC [23] | 0.733 | 0.076 | 0.790
OccNet [25] | 0.772 | 0.082 | 0.799
ConvOccNet [29] | 0.870 | 0.048 | 0.933
Ours | 0.875 | 0.043 | 0.945
Table 2: Quantitative comparison for point clouds reconstruction.

###  4.3 Voxel Super-Resolution

We also evaluate our method on the voxel super-resolution task, which takes
coarse voxelizations ($32^{3}$) as input, and reconstruct meshes with high-
resolution details. We use the same backbone as [29]. The quantitative results
are shown in Table 3. Our method achieves better performance compared with
other baselines.

It is noted that the coarse-to-fine expansion point strategy is simplified in
the voxel super-resolution task as the units near the surface can be obtained
from the voxel directly. For a unit of the input voxel, it is identified as a
unit near the surface if the number of units near it is not 26 or 0. We
observe that our approach can improve the performance comprehensively to other
baselines.

###  4.4 Inference Speed

To demonstrate the inference speed promotion of our approach, we compare with
OccNet [25] and ConvOccNet [29] on single-view and point cloud reconstruction.
For our approach, once the expansion points are sampled and the corresponding
Taylor series coefficients are predicted, the SDF value of arbitrary points
com be computed with the Taylor series of nearby expansion points. In
addition, the number of sampled expansion points from the coarse-to-fine
strategy will not be influenced by the resolution of mesh generation. So,
after the expansion points are obtained, the inference of our approach is
independent of the deep neural network. In this section, we aim to statistic
and analyze the inference speed on both the expansion sampling and the point
evaluation.

We compare the inference speed of our approach with OccNet [25] and ConvOccNet
[29] on mesh generation. The backbone of our approach on single-view and point
cloud data types are the same as OccNet and ConvOccNet respectively. All the
experiments in this section are implemented on one Titan V GPU. We only report
the time spent on evaluating points. The time spent on marching cubes is not
considered.

|  IoU $\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score $\uparrow$
---|---|---|---
Input | 0.631 | 0.136 | 0.440
OccNet [25] | 0.703 | 0.110 | 0.656
ConvOccNet [29] | 0.752 | 0.091 | 0.729
Ours | 0.761 | 0.081 | 0.755
Table 3: Quantitative comparison for voxel super-resolution. The first row is
the results of transforming the input voxels into meshes and computing the
metrics directly.

We can observe that, as the mesh resolution increases, the time increase of
our approach is less significant than other baselines. As the mesh resolution
increases, the time spent on other baselines increases about 4 to 8 times
while our approach only increases by about 2 times. This is because our
coarse-to-fine strategy for sampling expansion points is independent of the
mesh resolution, the time on sampling expansion points would not increase with
the mesh resolution increasing.

Besides, our approach spends more time on single-view reconstruction than
point cloud reconstruction. The reason is that identifying the units near the
surface is easier on point cloud data type. To guarantee the units which are
identified near the surface can cover the surface as much as possible, the
threshold $\epsilon_{in}$ and $\epsilon_{out}$ on single-view need to preserve
a larger margin than point cloud. Based on this, the number of expansion
points sampled in single-view reconstruction is often larger than point cloud
reconstruction.

Resolution | Single-view | Point cloud
---|---|---
OccNet | Ours | OccNet | ConvOccNet | Ours
64 | 0.101 | 0.031 | 0.171 | 0.113 | 0.024
128 | 0.783 | 0.069 | 0.916 | 0.492 | 0.035
Table 4: We statistic the average time spent (second) for each input sample in
the point evaluation step. Figure 6: Qualitative results on point cloud
reconstruction task with different Taylor series orders. Adopting Taylor
series of order 2 or 3 has achieved excellent reconstruction results compared
with the ground truth mesh. We do not visualize the mesh generated with order
4 since it has almost no difference from order 3.

###  4.5 Ablation Studies

In this section, we evaluate how the different orders of the Taylor series and
different numbers of nearest expansion points in the inference affect the
performance.

Order of Taylor Series To test the influence of the Taylor series with
different orders, we visualize the generated meshes which are evaluated with
coefficients of different orders. For example, when we evaluate the results of
the Taylor series with the order equal to 2, we only preserve the coefficients
of terms with the order less than or equal to 2. The coefficients of Taylor
series with the order larger than 2 are set to 0.

The qualitative results are shown in Figure 6. We observe that the terms with
the order 0 and 1 of the Taylor series can provide a rough outline of the
object while the terms with the higher order such as 2 and 3 represent the
fine details of the object. Besides, this result also supports the
aforementioned conclusion in Figure 1 (b), _i.e_., the low-order Taylor series
can provide an excellent shape approximation result. The improvement of
increasing the order of the Taylor series from order 3 is too slight to be
observed, which is also indicated by the quantitative results in Table 5.

Order |  IoU $\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score $\uparrow$
---|---|---|---
0 | 0.726 | 0.101 | 0.762
1 | 0.800 | 0.057 | 0.890
2 | 0.874 | 0.043 | 0.944
3 | 0.875 | 0.043 | 0.945
4 | 0.875 | 0.043 | 0.945
Table 5: The quantitative comparison of using different orders of Taylor
series to reconstruct from point clouds. Figure 7: The visualization of the
reconstruction results from a ”plane” point cloud with varying numbers of
expansion points at inference phase.

Number of Nearest Expansion Points in Inference As shown in Equation 5, we
compute the weighted average of the Taylor series with $k$-nearest expansion
points. To analyze how the value of $k$ influences the performance, we set
$k=1,2,4,8$ respectively in the inference phase, the results are shown in
Table 6.

From the results, we can see that the weighted average operation among
multiple expansion points promotes the performance significantly compared with
using a single expansion point when $k\leq 4$. However, when $k>4$, the
performance degrades slightly because the expansion points far away the target
query points provide less accurate prediction. Besides, exploiting larger $k$
will lead to more memory usage. We empirically find $k=4$ is a cost-effective
setting and adopt it in all experiments.

To better demonstrate the influence of exploiting $k$-nearest expansion
points, we visualize the reconstruction results of the ”plane” point cloud
with different values of $k$ in Figure 7. It can be observed that by using a
proper number of expansion points and performing the weighted average, we can
generate a smoother mesh and reduce the topology artifacts significantly.

k | Single-view | Point cloud
---|---|---
IoU$\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score$\uparrow$ |  IoU$\uparrow$ |  Chamfer-$L_{1}\downarrow$ |  F-Score$\uparrow$
1 | 0.593 | 0.193 | 0.543 | 0.867 | 0.045 | 0.940
2 | 0.598 | 0.192 | 0.551 | 0.872 | 0.044 | 0.944
4 | 0.599 | 0.192 | 0.553 | 0.875 | 0.043 | 0.945
8 | 0.600 | 0.196 | 0.550 | 0.876 | 0.043 | 0.946
Table 6: The influence of the number of nearest expansion points.

##  5 Conclusion

In this work, we propose TaylorImNet, a compact implicit 3D shape
representation based on the Taylor series for fast 3D shape reconstruction. In
contrast to classical deep implicit functions, our Taylor-based implicit
function is independent of the network in the inference phase, which
significantly boosts the mesh generation speed than other baselines.

Our experiments on different input data types demonstrate that our TaylorImNet
can not only speed up the mesh generation but also have comparable
representation capacity compared with the existing state-of-the-art baselines.
This proves the efficient representation capacity of our approach.

Broader Impact and Limitation Our experiments exploit the ShapeNet dataset
which is publicly available, thus do not involve privacy issues. However, we
did not test our approach on human datasets, thus cannot guarantee the
generalization ability to model humans. Besides, TaylorImNet requires sampling
a group of query points around each expansion point at the training phase,
which leads to a larger amount of training query points than the baselines
such as [25] and [29].

## References

  * [1] Panos Achlioptas, Olga Diamanti, Ioannis Mitliagkas, and Leonidas Guibas.  Learning representations and generative models for 3d point clouds.  In International conference on machine learning, pages 40–49. PMLR, 2018.
  * [2] Marc Alexa, Johannes Behr, Daniel Cohen-Or, Shachar Fleishman, David Levin, and Claudio T. Silva.  Computing and rendering point set surfaces.  IEEE Transactions on visualization and computer graphics, 9(1):3–15, 2003.
  * [3] Matan Atzmon and Yaron Lipman.  Sal: Sign agnostic learning of shapes from raw data.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [4] Rohan Chabra, Jan E Lenssen, Eddy Ilg, Tanner Schmidt, Julian Straub, Steven Lovegrove, and Richard Newcombe.  Deep local shapes: Learning local sdf priors for detailed 3d reconstruction.  In European Conference on Computer Vision, pages 608–625. Springer, 2020.
  * [5] Angel X Chang, Thomas Funkhouser, Leonidas Guibas, Pat Hanrahan, Qixing Huang, Zimo Li, Silvio Savarese, Manolis Savva, Shuran Song, Hao Su, et al.  Shapenet: An information-rich 3d model repository.  arXiv preprint arXiv:1512.03012, 2015.
  * [6] Zhiqin Chen and Hao Zhang.  Learning implicit fields for generative shape modeling.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5939–5948, 2019.
  * [7] Julian Chibane, Thiemo Alldieck, and Gerard Pons-Moll.  Implicit functions in feature space for 3d shape reconstruction and completion.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6970–6981, 2020.
  * [8] Christopher B Choy, Danfei Xu, JunYoung Gwak, Kevin Chen, and Silvio Savarese.  3d-r2n2: A unified approach for single and multi-view 3d object reconstruction.  In European conference on computer vision, pages 628–644. Springer, 2016.
  * [9] Yu Deng, Jiaolong Yang, and Xin Tong.  Deformed implicit field: Modeling 3d shapes with learned dense correspondence.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10286–10296, 2021.
  * [10] Yueqi Duan, Haidong Zhu, He Wang, Li Yi, Ram Nevatia, and Leonidas J Guibas.  Curriculum deepsdf.  In European Conference on Computer Vision, pages 51–67. Springer, 2020.
  * [11] Haoqiang Fan, Hao Su, and Leonidas J Guibas.  A point set generation network for 3d object reconstruction from a single image.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 605–613, 2017.
  * [12] Kyle Genova, Forrester Cole, Avneesh Sud, Aaron Sarna, and Thomas Funkhouser.  Local deep implicit functions for 3d shape.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4857–4866, 2020.
  * [13] Kyle Genova, Forrester Cole, Daniel Vlasic, Aaron Sarna, William T Freeman, and Thomas Funkhouser.  Learning shape templates with structured implicit functions.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 7154–7164, 2019.
  * [14] Rohit Girdhar, David F Fouhey, Mikel Rodriguez, and Abhinav Gupta.  Learning a predictable and generative vector representation for objects.  In European Conference on Computer Vision, pages 484–499. Springer, 2016.
  * [15] Amos Gropp, Lior Yariv, Niv Haim, Matan Atzmon, and Yaron Lipman.  Implicit geometric regularization for learning shapes.  arXiv preprint arXiv:2002.10099, 2020.
  * [16] Thibault Groueix, Matthew Fisher, Vladimir G Kim, Bryan C Russell, and Mathieu Aubry.  A papier-mâché approach to learning 3d surface generation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 216–224, 2018.
  * [17] Christian Häne, Shubham Tulsiani, and Jitendra Malik.  Hierarchical surface prediction for 3d object reconstruction.  In 2017 International Conference on 3D Vision (3DV), pages 412–420. IEEE, 2017.
  * [18] Zekun Hao, Hadar Averbuch-Elor, Noah Snavely, and Serge Belongie.  Dualsdf: Semantic shape manipulation using a two-level representation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7631–7641, 2020.
  * [19] Chiyu Jiang, Avneesh Sud, Ameesh Makadia, Jingwei Huang, Matthias Nießner, Thomas Funkhouser, et al.  Local implicit grid representations for 3d scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6001–6010, 2020.
  * [20] Michael Kazhdan, Matthew Bolitho, and Hugues Hoppe.  Poisson surface reconstruction.  In Proceedings of the fourth Eurographics symposium on Geometry processing, volume 7, 2006.
  * [21] Michael Kazhdan and Hugues Hoppe.  Screened poisson surface reconstruction.  ACM Transactions on Graphics (ToG), 32(3):1–13, 2013.
  * [22] Nikos Kolotouros, Georgios Pavlakos, and Kostas Daniilidis.  Convolutional mesh regression for single-image human shape reconstruction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4501–4510, 2019.
  * [23] Yiyi Liao, Simon Donne, and Andreas Geiger.  Deep marching cubes: Learning explicit surface representations.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2916–2925, 2018.
  * [24] William E Lorensen and Harvey E Cline.  Marching cubes: A high resolution 3d surface construction algorithm.  ACM siggraph computer graphics, 21(4):163–169, 1987.
  * [25] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger.  Occupancy networks: Learning 3d reconstruction in function space.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4460–4470, 2019.
  * [26] Jiteng Mu, Weichao Qiu, Adam Kortylewski, Alan Yuille, Nuno Vasconcelos, and Xiaolong Wang.  A-sdf: Learning disentangled signed distance functions for articulated shape representation.  arXiv preprint arXiv:2104.07645, 2021.
  * [27] Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove.  Deepsdf: Learning continuous signed distance functions for shape representation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 165–174, 2019.
  * [28] Despoina Paschalidou, Ali Osman Ulusoy, and Andreas Geiger.  Superquadrics revisited: Learning 3d shape parsing beyond cuboids.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [29] Songyou Peng, Michael Niemeyer, Lars Mescheder, Marc Pollefeys, and Andreas Geiger.  Convolutional occupancy networks.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part III 16, pages 523–540. Springer, 2020.
  * [30] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [31] Charles R Qi, Li Yi, Hao Su, and Leonidas J Guibas.  Pointnet++: Deep hierarchical feature learning on point sets in a metric space.  arXiv preprint arXiv:1706.02413, 2017.
  * [32] Daxuan Ren, Jianmin Zheng, Jianfei Cai, Jiatong Li, Haiyong Jiang, Zhongang Cai, Junzhe Zhang, Liang Pan, Mingyuan Zhang, Haiyu Zhao, et al.  Csg-stump: A learning friendly csg-like representation for interpretable shape parsing.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 12478–12487, 2021.
  * [33] Shunsuke Saito, Zeng Huang, Ryota Natsume, Shigeo Morishima, Angjoo Kanazawa, and Hao Li.  Pifu: Pixel-aligned implicit function for high-resolution clothed human digitization.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 2304–2314, 2019.
  * [34] Mo Shan, Qiaojun Feng, You-Yi Jau, and Nikolay Atanasov.  Ellipsdf: Joint object pose and shape optimization with a bi-level ellipsoid and signed distance function description.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5946–5955, 2021.
  * [35] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Advances in Neural Information Processing Systems, 33, 2020.
  * [36] Matthew Tancik, Pratul P Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  arXiv preprint arXiv:2006.10739, 2020.
  * [37] Nanyang Wang, Yinda Zhang, Zhuwen Li, Yanwei Fu, Wei Liu, and Yu-Gang Jiang.  Pixel2mesh: Generating 3d mesh models from single rgb images.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 52–67, 2018.
  * [38] Jiajun Wu, Chengkai Zhang, Tianfan Xue, William T Freeman, and Joshua B Tenenbaum.  Learning a probabilistic latent space of object shapes via 3d generative-adversarial modeling.  arXiv preprint arXiv:1610.07584, 2016.
  * [39] Haozhe Xie, Hongxun Yao, Xiaoshuai Sun, Shangchen Zhou, and Shengping Zhang.  Pix2vox: Context-aware 3d reconstruction from single and multi-view images.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 2690–2698, 2019.
  * [40] Shun Yao, Fei Yang, Yongmei Cheng, and Mikhail G Mozerov.  3d shapes local geometry codes learning with sdf.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 2110–2117, 2021.
  * [41] Kangxue Yin, Zhiqin Chen, Hui Huang, Daniel Cohen-Or, and Hao Zhang.  Logan: Unpaired shape transform in latent overcomplete space.  ACM Transactions on Graphics (TOG), 38(6):1–13, 2019.
  * [42] Zerong Zheng, Tao Yu, Qionghai Dai, and Yebin Liu.  Deep implicit templates for 3d shape representation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1429–1439, 2021.


# Multiphoton Spectroscopy of a Dynamical Axion Insulator

Olivia Liebman<EMAIL_ADDRESS>Jonathan Curtis  College of Letters and
Science, University of California, Los Angeles  Ioannis Petrides  College of
Letters and Science, University of California, Los Angeles  Prineha Narang
<EMAIL_ADDRESS>College of Letters and Science, University of California, Los
Angeles

###### Abstract

The unusual magnetoelectric transport present in Weyl semimetals can be
compactly understood as manifestations of an underlying axion field, which
itself is determined by the microscopic band structure. The axion couples
nonlinearly to electric and magnetic fields and possesses a signature
topological magnetoelectric response, leading to a modified form of Maxwell’s
equations known as axion electrodynamics. Axions are naturally present in Weyl
semimetals relating to the separation of the Weyl nodes in energy and in
crystal momentum. In the presence of strong interactions, charge density-wave
(CDW) order may develop which serves to gap the Weyl nodes and introduces
corresponding collective excitations. When the inherent chiral symmetry of
Weyl semimetals is spontaneously broken by the formation of CDW order, the
resultant chiral condensate is endowed with intrinsic dynamics which yields a
dynamical contribution to the axion field. However, unambiguous identification
of this dynamical axion mode is challenging due to its inherent nonlinear
coupling to electromagnetic fields. Therefore, we propose an all-optical
protocol for verifying and characterizing dynamical axion collective modes in
Weyl semimetals with CDW order. First, we show that axion collective mode can
be excited using two-photon excitation schemes. Following excitation, the
collective axion oscillations are then diagnosed by measuring the time-
resolved Kerr rotation. Our results demonstrate a pathway towards utilizing
multi-photon and entangled pair spectroscopies to identify new correlated
phases in quantum matter.

††preprint: APS/123-QED

Introduction—Axionic particles were originally proposed in high-energy physics
in order to solve the charge-parity problem as the Nambu-Goldstone bosons
associated with a new global axial U(1) symmetry [1, 2]. Today, axions are a
leading dark matter candidate and could help explain the matter-antimatter
asymmetry present in the universe [3]. While the axion has yet to be
discovered in particle physics, its condensed matter [4] counterpart has been
theorized to exist in certain three-dimensional materials [5, 6, 7, 8, 9]. A
characteristic identifier of the presence of an axion is a signature
magnetoelectric response, leading to a modified form of electrodynamics known
as axion electrodynamics [10, 11, 12].

It turns out Weyl semimetals naturally have an axion term that is related to
the separation in k-space and in energy of the Weyl nodes. What’s more, when
the Weyl system undergoes lattice translation symmetry breaking due to the
emergence of charge density wave (CDW) order [13, 14], the collective motion
associated to the phase of the CDW in turn leads to a dynamical axion effect
and to a nonlinear magnetoelectric effect on top of the response due to band
structure effects [15, 16, 17, 18]. In the presence of uncompensated carrier
densities this may even lead to spatially inhomogeneous textures due to
softening of the axionic collective mode [19].

In addition to being of fundamental interest, there have also been proposals
that the magnetoelectric coupling mediated by the axion response in materials
may be of practical value, e.g., in nonreciprocal thermal emitters or
rectifiers [20, 21]. However, while some progress has been made, techniques
for the unambiguous identification of dynamical axion quasiparticles are
challenging, and often rely on studying linear perturbations in the presence
of a strong magnetic field [22, 5, 23, 24].

Here we lay out an all-optical, contactless protocol for identification of
nonlinear signatures of dynamical axions in a Weyl semimetal with CDW order.
This protocol breaks down into two parts; first a dynamical axion mode is
nonlinearly excited through two-photon processes, which carries a distinct
dependence on the incident angle, polarization, and frequencies of the two
beams. Following excitation, the induced axion dynamics are then detected
through their manifestation in the Kerr angle rotation of a third reflected
probe beam, which acts as a faithful proxy for the dynamical evolution of the
axion in the material [25, 26]. More generally, our work demonstrates the
potential of using entanglement pairs to identify correlated phases in quantum
materials.

Axion electrodynamics—We now describe the excitation protocol, starting from a
model for dynamical axion collective modes in a Weyl semimetal with CDW order.
The simplest form to discuss Weyl physics is when time-reversal symmetry is
broken but inversion symmetry is persevered, thus there are minimally two
nodes of equal energy. The CDW Weyl Hamiltonian expanded near the degeneracy
points is

$\displaystyle\bm{\mathcal{H}}(\bm{k})=\left(\begin{array}[]{cc}v_{F}\bm{\sigma}\cdot\bm{k}&\Delta\\\
\bar{\Delta}&-v_{F}\bm{\sigma}\cdot\bm{k}\end{array}\right)$ (3)

where $v_{F}$ is the Fermi velocity, $\bm{\sigma}$ is the triplet of Pauli
matrices, $\bm{k}$ is the crystal momentum, and $\Delta=|\Delta|e^{i\theta}$
is the CDW term and is analogous to a complex mass term, also represented in
terms of the gap amplitude $|\Delta|$ and phase $\theta$. In an uncorrelated
Weyl semimetal with no CDW $\Delta=0$ the nodes are gapless. The chirality of
the Weyl nodes is identical to the topological charge associated to each band
crossing, given by the Chern flux $C=+1/2$ ($C=-1/2$) for the right- (left-)
handed Weyl particles [27]. This chiral charge can be understood as magnetic
monopoles, or sources and sinks, of the Berry curvature which is singular at
the nodes and acts as an effective magnetic field in reciprocal space [28].

In the presence of sufficiently strong interactions, a transition to CDW
ordered phase with complex order parameter $\Delta$ can occur. If the order is
incommensurate, by Goldstone’s theorem this will have a soft sliding mode
associated with the phase of the CDW order, $\theta(x)$. Due to the chiral
anomaly, it is known that the phase of the complex CDW can be identified with
the axion response [6] of the material. This serves to modify the
electrodynamics of the medium, with Lagrangian

$\displaystyle\mathcal{L}_{\rm
EM}=\frac{1}{2}(\epsilon\bm{E}^{2}-\frac{1}{\mu}\bm{B}^{2})+g\theta\bm{E}\cdot\bm{B}.$
(4)

The first two terms are the sourceless Maxwell Lagrangian in terms of the
electric (magnetic) fields $\bm{E}$ ($\bm{B}$), and $\epsilon$ ($\mu$) is the
static dielectric (permeability) tensor (we use units with
$\epsilon_{0}=\mu_{0}=1)$. The last term proportional to $\bm{E}\cdot\bm{B}$
is the so-called axion response that leads to a magnetoelectric coupling
mediated by the CDW phase $\theta$, with coupling constant $g=\alpha/\pi$.

The total CDW phase $\theta$ is given by

$\displaystyle\theta(\bm{r},t)=\bm{Q}\cdot\bm{r}+\delta\theta(\bm{r},t)\,,$
(5)

where the first term is dubbed the static axion field as it depends only on
the separation in momentum-space $\bm{Q}$ between the Weyl nodes, and is
determined by the band-structure of the material 111In some cases it is also
possible to have a separation of the Weyl nodes in energy of $Q_{0}$, leading
to an additional contribution of $Q_{0}t$ in the static axion. This leads to
additional topological responses but we will not consider this case here.. The
second term $\delta\theta(\bm{r},t)$ is the dynamical axion field, which can
be understood as the collective mode associated to the sliding phase mode of
the CDW.

The electromagnetic Lagrangian (4) is supplemented by the Nambu-Goldstone
equations of motion of the collective dynamical axion mode, namely

$\displaystyle\mathcal{L}_{\rm
NG}=\frac{\kappa}{2}[(\partial_{t}\delta\theta)^{2}-v^{2}(\bm{\nabla}\delta\theta)^{2}-\Omega_{0}^{2}\delta\theta^{2}]$
(6)

The first term in the above equation is the kinetic energy with $\kappa$ the
chiral compressibility, which is proportional to the density of states at the
Fermi level; the second term describes the spatial dispersion with $v$ the
speed of sound in the medium; and the final term is the pinning-induced gap
term with $\Omega_{0}$ the pinning frequency.

Figure 1:  Schematic of the general excitation and detection set-up. Two
incident beams come in at a relative angle of $\pi-2\alpha$ with polarizations
$\bm{\hat{\varepsilon}}_{1},\bm{\hat{\varepsilon}}_{2}$; frequencies
$\omega_{1},\omega_{2}$; and wave vectors $\bm{q}_{1},\bm{q}_{2}$; and are
used to induce finite $\bm{E}\cdot\bm{B}$, exciting an axionic collective mode
using a multiphoton absorption event. In order to verify the creation of the
axionic mode, time-resolved Kerr rotation can be used. For a probe beam at
normal incidence (assumed to be collinear with the static axion response along
$\bm{Q}$), the reflected beam’s polarization will be elliptically polarized at
an angle that oscillates in phase with the induced axion oscillations. Figure
2:  Induced axion amplitude $|\delta\theta|^{2}$ as a function of incident
photon frequencies $\omega_{1},\omega_{2}$ for different incident angles
$\alpha$ in the maximally-crossed polarization channel. (a) Excitation
amplitude for $\alpha=0$ (beams anti-parallel head on). We see a pronounced
two-photon absorption feature when $\omega_{2}+\omega_{1}=\Omega_{0}$. (b)
Orthogonal incidence when $\alpha=\pi/4$. We see that in addition to the two-
photon absorption channel, we also see the emergence of a much stronger
inelastic (stimulated Raman) excitation channel. Note the $y=0$ line separates
the two half planes, as well as note the separate scale bars for each. This
was done to resolve the axion response in the upper half plane, which is
significantly weaker as compared to the response shown in the lower half
plane. (c) Parallel collinear incidence ($\alpha=\pi/2$). In this case, the
two-photon absorption channel has been completely suppressed and instead the
Raman excitation channel is maximal. Note in this case the relative momentum
of the two beams is zero and therefore there must be a finite relative-
frequency.

The equations of motion of the gauge potential derived from Eq. (6) are given
by

$\partial^{2}_{t}\bm{A}-\bm{\nabla}\times\bm{\nabla}\times\bm{A}=g(\partial_{t}\theta\bm{\nabla}\times\bm{A}-\bm{\nabla}\theta\times\partial_{t}\bm{A})$
(7)

where $\bm{B}=\bm{\nabla}\times\bm{A}$ and $\bm{E}=\partial_{t}\bm{A}$ are the
magnetic and electric fields, respectively, given in terms of the gauge
potential with the Weyl gauge selected $\bm{\nabla}\phi=0$ for scalar
potential $\phi$. The current response of the axion insulator is hence

$\bm{J}=g\bm{\nabla}\theta\times\bm{E}+g\dot{\theta}\bm{B}.$ (8)

In particular, the first term is equivalent to an anomalous Hall response
proportional to the spatial gradient of the axion field, while the second term
is known as the chiral magnetic effect and is only relevant in the presence of
time-dependent axion fields. The equations of motion for the dynamical axion
$\delta\theta$ are

$\kappa(\partial_{t}^{2}+\gamma\partial_{t}-v^{2}\nabla^{2}+\Omega_{0}^{2})\delta\theta=g\bm{E}\cdot\bm{B}$
(9)

Here the second term introduces a phenomenological damping parameter $\gamma$.
These equations qualitatively describe massive sound waves driven by an
effective force $g\bm{E}\cdot\bm{B}$. Thus as long as the electric and
magnetic fields have some non-orthogonal component the collective dynamical
axion field can be excited.

Exciting the dynamical axion mode—As seen in Eq. (9), exciting a dynamical
axion mode in a material is a nonlinear optical process requiring at least two
electromagnetic fields. In general, the intensity of the axion response
depends on the frequency, momenta, polarization, and angle of the incident
beams, see Fig. 1 for proposed protocol. We consider the following plane-wave
ansatz $\bm{A}(q)=\bm{A}_{1}\delta_{q,q_{1}}+\bm{A}_{2}\delta_{q,q_{2}}$,
where $q_{i}=(\omega_{i},\bm{q}_{i})$ is the four-momentum with $\omega_{i}$
the frequency and $\bm{q}_{i}$ the momentum. The field amplitudes
$\bm{A}_{i}=A_{i}\bm{\hat{\varepsilon}}_{i}$ are related to the polarizations
$\bm{\hat{\varepsilon}}_{i}$ which obey
$\bm{q}_{i}\cdot\bm{\hat{\varepsilon}}_{i}=0$. We further assume that the
momenta satisfy the free-space dispersion relations so that
$|\bm{q}_{i}|=|\omega_{i}|$.

The response of the dynamical axion to the application of the vector potential
$\bm{A}(q)$ is given by

$\displaystyle\delta\theta(\omega_{1},\omega_{2},\bm{q}_{1},\bm{q}_{2})=$
$\displaystyle\frac{(g/\kappa)(\omega_{2}\bm{q}_{1}-\omega_{1}\bm{q}_{2})\cdot(\bm{\hat{\varepsilon}_{1}}\times\bm{\hat{\varepsilon}_{2})}}{(\omega_{1}+\omega_{2})^{2}+i(\omega_{1}+\omega_{2})\gamma-\Omega_{0}^{2}}A_{1}A_{2}$
(10)

here the term proportional to $v^{2}\nabla^{2}$ of Eq. (9) has been dropped
since $v\ll c$. Importantly, the appearance of the vector triple product
constrains the amplitude of the dynamical axion mode to the geometric volume
defined by the two polarization vectors $\hat{\bm{\varepsilon}}_{i}$ and the
relative four-momenta $q_{i}$ of the incident electromagnetic fields; as a
result, this provides a route towards smoking-gun identification of the axion
collective mode.

In order to excite the axion, the polarizations must be non-collinear and also
must be non-collinear with the relative four-momentum
$\omega_{1}\bm{q}_{2}-\omega_{2}\bm{q}_{1}$. In Fig. 2 we explore the efficacy
of this two-photon excitation protocol for a variety of different incident
beam angles $\alpha$ and frequencies $\omega_{i}$ in the maximally-crossed
polarization channel, c.f., Fig. 1. Depending on the relative angle $\alpha$
there are two two-photon channels which can be used to excite the dynamical
axion mode. Two-photon absorption can be used when the beams are close to
head-on, i.e., $\alpha=0$. For parallel-propagating collinear beams, with
$\alpha=\pi/2$ this channel is closed, and instead stimulated Raman excitation
can be used to excite the dynamical axion mode by inelastic scattering of
light. At intermediate incident angles, e.g., $\alpha=\pi/4$, both channels
are active. Fig. 2 shows clearly how varying degrees of freedom of the
incident beam will induce changes to axion excitation intensity, with the
inelastic stimulated Raman channel for collinear beams resulting in a stronger
response relative to the two-photon absorption process.

Axion detection via Kerr angle modulation—Now that the protocol for exciting
an axion collective mode has been established, we move on demonstrate how it
can be detected via an optical signature. The key idea is to use a linearly
polarized detection beam of known frequency at normal incidence, to allow for
optical detection of the axion field by its effect on the reflected light; a
technique well established for the spectroscopy of static axion responses
[25].

The proposed pump-probe style protocol is illustrated in Fig. 1, where the
pump beam excites the dynamical axion mode in the bulk and is bounded above by
a vacuum with relative permittivity $\epsilon=1$. For simplicity we assume the
material is nonmagnetic so that $\mu=\mu_{0}=1$. Also, by assuming the probe
incident beam wave vector along the z-axis, it is sufficient to consider the
bulk dielectric permittivity tensor only in the plane orthogonal to the
direction of propagation of the wave

$\epsilon(\omega)=\begin{pmatrix}\epsilon_{xx}&\epsilon_{xy}\\\
\epsilon_{yx}&\epsilon_{yy}\end{pmatrix}=\begin{pmatrix}\epsilon_{1}&i\epsilon_{2}\\\
-i\epsilon_{2}&\epsilon_{1}\end{pmatrix}$ (11)

where the dielectric tensor includes Hall effect contributions given as the
complex off-diagonal matrix elements. Furthermore, we assume that the axion
field’s spatial variation is small compared to the wavelength of the incident
light, such that $\bm{\nabla}\delta\theta\ll\bm{q}\delta\theta$ everywhere
except close to the surface where we treat $\nabla\delta\theta$ as a sharp
$\delta$-function singularity. We also take the adiabatic approximation where
$\partial_{t}\delta\theta\ll\omega\delta\theta$ such that contributions from
the time derivative are sufficiently small so as to be taken to be zero, with
the justification being that the axion frequency varies slowly as compared to
the frequency of light. This will be later relaxed to allow for dynamical
axion modes.

Figure 3:  Manifestations of excited axion collective mode in transient Kerr
angle rotation. Using protocol established here the axion is pumped by a
finite $\bm{E}\cdot\bm{{B}}$ profile, which is tuned to be in two-photon
resonance with the dynamical axion. This induces collective oscillations of
the axion angle shown above as $\delta\theta(t)$. These oscillations then
manifest in the time-dependence of the Kerr-angle $\Theta_{K}(t)$ which
approximately track the axion field. Parameters used here are: axion resonant
frequency $\Omega_{0}=1$, axion damping $\gamma=0.2$, constant $g=0.1$, chiral
compressiblility $\kappa=1$, probe frequency $\omega=1$, and dielectric matrix
elements $\epsilon_{1}=10$ and $\epsilon_{2}=0.7$. Note the amplitudes for
each response are not to scale. The pump amplitude is an order of magnitude
greater than the $\delta\theta$ amplitude, which is itself an order of
magnitude greater than the $\Theta_{K}$ amplitude.

Under these assumptions Eq. (7) reduces to
$\bm{\nabla}\times\bm{\nabla}\times\bm{A}=-\epsilon(\bm{r})\partial_{t}^{2}\bm{A}$
in the medium. Solving for the gauge potential $\bm{A}$ naturally yields the
chiral basis set $\bm{A}^{\pm}_{j}=A_{j}^{\pm}\frac{1}{\sqrt{2}}\left(1,\pm
i\right)^{T}$ where $A_{j}^{\pm}$ is an amplitude prefactor, and corresponding
eigenvalues $q^{\pm}=\omega n^{\pm}$ with
$n_{\pm}=\sqrt{\epsilon_{1}\pm\epsilon_{2}}$. The relfectance and
transmittance amplitudes can be obtained via the continuity of electric and
magnetic fields at the interface, given as $\bm{E}_{i}+\bm{E}_{r}=\bm{E}_{t}$,
and $\bm{B}_{i}+\bm{B}_{r}=\bm{B}_{t}+g\delta\theta\bm{E}_{t}$. We are most
interested in nonlinear optics effects which may be detected using
spectroscopic techniques. In turn, the coefficients of most use are the
reflection coefficients, given in terms of the incident beam as

$R^{\pm}=\pm\frac{1\mp(n^{\pm}\mp ig\delta\theta)}{1\pm(n^{\pm}\mp
ig\delta\theta)}I^{\pm}$ (12)

Here $I^{\pm}$ denotes the incident light written in the chiral basis. The
plus/minus reflection coefficients are, in general, different for left- and
right-circularly polarized light.

To promote $\delta\theta$ to time-dependence, the partial differential
equation of Eq. (9) are solved with the initial condition $\delta\theta(t\ll
t_{pump})=0$ as the axion field is not excited before the onset of the pump at
time $t_{pump}$. Following excitation, the damping term $\gamma$ causes a ring
down of the axion field, see Fig. 3. The key result is the reflection
amplitude’s direct dependence on the axion field. This will manifest as a
modulation of the reflected beam Kerr angle, and serves as our smoking gun
optical signature to unambiguously deduce the presence of a dynamical axion
field in the material.

We now do an order-of-magnitude estimate of the size of this effect. The
change in the polarization rotation depends on the size of the induced
$g\delta\theta$ (recall $g=\alpha/\pi$ is a universal constant). Applying the
equations of motion, we find that this will scale as
$g\delta\theta\sim(g^{2}/\kappa)(E/\Omega)^{2}/c$, as $B\sim E/c$ in these
units, and $\Omega$ is the dynamical axion resonance frequency which we take
to be 10GHz. Up to numerical factors we take $\kappa=\nu(E_{F})\sim
k_{F}^{2}/v_{F}$ for a semimetal with Fermi velocity $v_{F}\sim 10^{-3}c$.
$k_{F}$ in turn is set by the carrier density, which we take to be
$k_{F}\sim.1\AA^{-1}$. For a large but physically realizable electric field of
order 1MV/cm this gives a response of $g\delta\theta\sim
g^{2}v_{F}/cE^{2}/(\Omega k_{F})^{2}\sim 1.2$, which is a small but observable
response.

Conclusion—Beginning with a Weyl semimetal with charge density wave order, we
have shown how the separation of the Weyl nodes in tandem with the CDW-induced
phase fluctuations can give rise to a dynamical axion field. To model the
dynamical axion mode and determine a novel optical signature we solved the
nonlinear partial differential equations describing its dynamics and
established a two-step, pump-probe style protocol which excites this mode with
a two-photon excitation scheme, and subsequently detects it with a third probe
beam. The intensity of the axion excitation is stronger for parallel collinear
beams relative to other beam configurations, as can be understood from Eq.
(23) and is shown in Fig. 2(c). Our key result is reflected probe beam’s Kerr
angle is modulated such that it approximately tracks the axion field
oscillation.

A motivation of this work is to provide theoretical predictions that can drive
future experimental work, where potential material candidates could include
axion insulators MnBi2Te4 [30] and (TaSe4)2I [22, 31, 32]. This paper has
demonstrated the potential for mutlidimensional [33, 34, 35] and multiphoton
spectroscopy in the study of correlated topological materials and their
dynamics [36]. The Quantum Chromodynamics axion, possibly within reach of
experimental verification by current experiments, represents one of the most
compelling pathways beyond the Standard Model, by providing at the same time
an elegant solution to the strong CP problem and an excellent dark matter
candidate [3, 37]. Meanwhile, condensed matter axions arising in novel,
topological materials has lead to a more accessible playground for studying
the exotic phenomena associated with these particles that may provide insight
when mapped back to high energy physics and impact experimental prospects of
an axion discovery via axion-mediated forces. Our work also has implications
for condensed matter axion systems that are being used to detect their high
energy counterparts in direct detection dark matter experiments [38, 39].

###### Acknowledgements.

The authors would like to acknowledge fruitful discussions with Aaron Chou,
Fahad Mahmood, Soyeun Kim, Ankit Disa, and Yikun Wang. This work is entirely
supported by the Quantum Science Center (QSC), a National Quantum Information
Science Research Center of the U.S. Department of Energy (DOE). P.N.
gratefully acknowledges support from the Gordon and Betty Moore Foundation
grant number #8048 and from the John Simon Guggenheim Memorial Foundation
(Guggenheim Fellowship).

##  Appendix A Derivation of axion response function

Here we provide more detailed derivations of the results presented in the main
work. We begin by outlining the steps to excite the axion mode via two
externally applied light beams, and arrive at the axion response function. To
start, we show how the gauge field and axion equations of motion may be
derived from the full axion modified electromagnetic Lagrangian

$\displaystyle\mathcal{L}=\frac{1}{2}(\epsilon\bm{E}^{2}-\frac{1}{\mu}\bm{B}^{2})+g\theta\bm{E}\cdot\bm{B}+$
$\displaystyle\frac{\kappa}{2}[(\partial_{t}\delta\theta)^{2}-v^{2}(\nabla\delta\theta)^{2}-\Omega_{0}^{2}\delta\theta^{2}]$
(13)

by utilizing the Euler-Lagrange formalism. For $\bm{A}$ this is

$\frac{\partial}{\partial t}\left(\frac{\partial
L}{\partial(\partial_{t}\bm{A})}\right)+\nabla\cdot\left(\frac{\partial
L}{\partial(\nabla\bm{A})}\right)-\frac{\partial L}{\partial\bm{A}}=0$ (14)

Varying with respect to $\bm{A}$ yields its equation of motion

$\partial^{2}_{t}\bm{A}-\nabla\times\nabla\times\bm{A}=g(\partial_{t}\delta\theta\nabla\times\bm{A}-\nabla\delta\theta\times\partial_{t}\bm{A})$
(15)

Similarly for the axion field

$\frac{\partial}{\partial t}\left(\frac{\partial
L}{\partial(\partial_{t}\theta)}\right)+\nabla\cdot\left(\frac{\partial
L}{\partial(\nabla\theta)}\right)-\frac{\partial L}{\partial\theta}=0$ (16)

and the corresponding $\theta$ equation of motion

$\kappa(\partial_{t}^{2}+\gamma\partial_{t}-v^{2}\nabla^{2}+\Omega_{0}^{2})\delta\theta=g[\partial_{t}\bm{A}(\nabla\times\bm{A})]$
(17)

As stated in the main text, $v^{2}\nabla^{2}\theta$ may be dropped by
observing $v\ll c$. The goal is to solve for the axion response to an
externally applied electromagnetic field. This is first carried out by
perturbatively expanding for $\theta$ and $\bm{A}$ with small parameter
$\varepsilon$

$\displaystyle\bm{A}(\bm{x},t)=\epsilon\bm{A}_{1}+\epsilon^{2}\bm{A}_{2}+\epsilon^{3}\bm{A}_{3}+...$
(18)
$\displaystyle\theta(\bm{x},t)=\epsilon\theta_{1}+\epsilon^{2}\theta_{2}+\epsilon^{3}\theta_{3}+...$
(19)

Performing this expansion and plugging into the gauge field and axion
equations of motion, respectively, determines the axion field must be
minimally excited by two gauge fields. With this in mind, we solve for
$\delta\theta$ by Fourier transforming Eq. (17). Note the vector potential
Fourier transforms as $\bm{A}(x_{i})=\int_{q}e^{iq_{i}\cdot x}\bm{A}(q_{i})$;
and we denote
$\int_{q}=\int\frac{d^{3}q}{(2\pi)^{3}}\int\frac{d\Omega}{2\pi}$, and $q\cdot
x=\bm{q}\cdot\bm{x}-\Omega t$ in the exponential.

$\displaystyle\kappa\int d^{4}xe^{-iq\cdot
x}(\partial_{t}^{2}+\gamma\partial_{t}+\Omega_{0}^{2})\delta\theta(x)=$
$\displaystyle\kappa(-\Omega^{2}+i\Omega+\Omega_{0}^{2})\int d^{4}xe^{-iq\cdot
x}\delta\theta(x)=$
$\displaystyle\kappa(-\Omega^{2}+i\Omega+\Omega_{0}^{2})\delta\theta(q)$ (20)

Inverting the matrix $(-\Omega^{2}+i\Omega+\Omega_{0}^{2})\mathbb{I}$ to
isolate $\delta\theta(q)$ on the left-hand side:

$\displaystyle\delta\theta(q,\omega)=[-\Omega^{2}+i\Omega+\Omega_{0}^{2}]^{-1}(\frac{g}{\kappa})\times$
$\displaystyle\int d^{4}xe^{-iq\cdot
x}[\partial_{t}\bm{A}(x_{1})\cdot\nabla\times\bm{A}(x_{1})]$
$\displaystyle=[-\Omega^{2}+i\Omega+\Omega_{0}^{2}]^{-1}(\frac{g}{\kappa})\times$
$\displaystyle\int d^{4}xe^{-iq\cdot x}[\partial_{t}\int e^{iq_{1}\cdot
x}\bm{A}_{1}(q_{1})\cdot\nabla\times\int e^{iq_{2}\cdot x}\bm{A}_{1}(q_{2})]$
$\displaystyle=[-\Omega^{2}+i\Omega+\Omega_{0}^{2}]^{-1}(\frac{g}{\kappa})\int_{q_{1},q_{2}}\int_{q}d^{4}x\times$
$\displaystyle e^{-i(q\cdot x-q_{1}\cdot x-q_{2}\cdot
x)}(-i\Omega_{1}\bm{A}_{1}(q_{1}))(i\bm{q_{2}}\times\bm{A}_{1}(q_{2}))$ (21)

Making use of the exponential representation as the Dirac delta function $\int
d^{4}xe^{-i(q\cdot x-q_{1}\cdot x-q_{2}\cdot x)}=\delta_{q_{1}+q_{2},q}$
simplifies this to

$\displaystyle\theta_{2}(q)=[-\Omega^{2}+\Omega_{0}^{2}]^{-1}(\frac{g}{\kappa})\times$
$\displaystyle\int_{q_{1},q_{2}}\delta_{q_{1}+q_{2},q}\Omega_{1}\bm{A}_{1}(q_{1})\cdot(\bm{q_{2}}\times\bm{A}_{1}(q_{2})).$
(22)

Recall the field is minimally composed of two beams whose frequencies must add
up to axion resonance $\Omega=\omega_{1}+\omega_{2}$ in order to excite the
mode. Utilizing a plane-wave ansatz
$\bm{A}(q)=\bm{A}_{1}\delta_{q,q_{1}}+\bm{A}_{2}\delta_{q,q_{2}}$, with four-
momenta $q_{i}=(\omega_{i},\bm{q}_{i})$ for frequency $\omega_{i}$ and
momentum $\bm{q}_{i}$; the solution for the axion field is

$\displaystyle\delta\theta(\omega_{1},\omega_{2},\bm{q}_{1},\bm{q}_{2})=$
$\displaystyle\frac{(g/\kappa)(\omega_{2}\bm{q}_{1}-\omega_{1}\bm{q}_{2})\cdot(\bm{\hat{E}_{1}}\times\bm{\hat{E}_{2})}}{(\omega_{1}+\omega_{2})^{2}+i(\omega_{1}+\omega_{2})\gamma-\Omega_{0}^{2}}A_{1}A_{2}$
(23)

##  Appendix B Axion modified electrodynamics

Next consider the interface of a dynamical axion insulator (DAI) and a vacuum.
The side of the vacuum is governed by the sourceless Maxwell’s equations;
while the other side is the DAI governed by axion electrodynamics with
modified constitutive equations

$\displaystyle\bm{D}=\epsilon\bm{E}-g\theta\bm{B}$ (24)
$\displaystyle\bm{H}=\frac{1}{\mu}\bm{B}+g\theta\bm{E}$ (25)

where $\epsilon=\epsilon(\bm{r},t)$ is the complex dielectric tensor and $\mu$
= $\mu(\bm{r},t)$ is the magnetic permeability tensor which we set to 1 for
simplicity, as is the case in nonmagnetic materials. In the medium, modified
Ampere’s law gives

$\displaystyle\nabla\times\bm{H}=\bm{J}+\partial_{t}\bm{D}\rightarrow$ (26)
$\displaystyle\nabla\times(\frac{1}{\mu_{0}}\bm{B}+g\theta\bm{E})=\bm{J}+\partial_{t}(\epsilon\bm{E}-g\theta\bm{B})$
(27)

where $\bm{J}=\sigma\bm{E}$ and we include in the conductivity $\sigma$ the
contribution from the Hall conductivity. Substituting back in for the vector
potential gives and simplifying
$\frac{1}{\mu}\nabla\times\nabla\times\bm{A}=\epsilon(\bm{r})\partial_{t}^{2}\bm{A}$
and so far it is assumed $\theta$ does not spatially vary. As stated in the
main text, we take the adiabatic approximation for now and let
$\dot{\theta}=0$. Later $\theta$ will be promoted to be time-dependent.
Fourier-transforming equation this expression

$\displaystyle\frac{1}{\mu}(i\bm{k}\times
i\bm{k}\times\bm{A})=\epsilon(\omega)\omega^{2}\bm{A}\rightarrow$ (28)
$\displaystyle\frac{1}{\mu}(\bm{k}\cdot\bm{k}\mathbb{I}-\bm{k}\otimes\bm{k})\bm{A}=\epsilon(\omega)\omega^{2}\bm{A}$
(29)

To simplify this slightly, let the incident and reflected light propagate in
the $\hat{z}$-direction and let this be orthogonal to the interface whose
surface lies in the $xy$-plane. Then due to boundary conditions the wave-
vector of the transmitted light will also be along the $\hat{z}$-direction so
that $\bm{k}=(0,0,k_{z})$. This simplifies the above matrix equation as

$\begin{pmatrix}k_{z}\omega^{2}&0\\\
0&k_{z}\omega^{2}\end{pmatrix}\begin{pmatrix}A_{x}\\\
A_{y}\end{pmatrix}=\omega^{2}\begin{pmatrix}\epsilon_{xx}&\epsilon_{xy}\\\
\epsilon_{yx}&\epsilon_{yy}\end{pmatrix}\begin{pmatrix}A_{x}\\\
A_{y}\end{pmatrix}$ (30)

Let the permittivity tensor have the form for a gyroelectric medium at
frequency $\omega$

$\epsilon(\omega)=\begin{pmatrix}\epsilon_{1}(\omega)&i\epsilon_{2}(\omega)&0\\\
-i\epsilon_{2}(\omega)&\epsilon_{1}(\omega)&0\\\
0&0&\epsilon_{3}(\omega)\end{pmatrix}$ (31)

Where we can ignore the z-component $\epsilon_{3}(\omega)$ which is zero for a
plane-wave propagating in the z-direction. Solving the matrix equation 30 the
eigenvalues correspond to
$\bm{k}^{\pm}=\omega\sqrt{\epsilon_{1}\pm\epsilon_{2}}\hat{e}^{\pm}$ and the
eigenvectors form a chiral basis
$\bm{A}_{j}^{\pm}=A_{j}^{\pm}\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
{\pm}i\end{pmatrix}$. The incident field may be rewritten in the new chiral
basis

$\displaystyle\bm{A}_{i}=\bm{A}_{i}^{+}\hat{e}_{+}+\bm{A}_{i}^{-}\hat{e}_{-}$
(32)

similarly for reflected and transmitted light.

Now let’s promote $\theta$ to be time-dependent. Solving the nonlinear partial
differential equation given in Eq. (17) for time varying electromagnetic
fields yields a dynamical contribution $\delta\theta$ to the axion field. This
requires an expression for the pump beam $p(t)\propto\bm{E}\cdot\bm{B}$ to
drive the axion mode, taken as a Gaussian wavepacket
$p(t)=e^{(t/\sigma)^{2}}cos(\Omega_{0}t)$, where $\sigma$ tunes the spread of
the Gaussian, and $\Omega_{0}$ denotes the axion resonant frequency. Thus the
axion term will be a dynamical fluctuating field, which we later use to
demonstrate how such a dynamical axion field generates a time-dependent
modulation of the Kerr rotation.

##  Appendix C Boundary conditions and Kerr rotation angle

The boundary conditions for electric and magnetic fields at an interface with
no surface charge or current are

$\displaystyle\bm{E}_{i}+\bm{E}_{r}=\bm{E}_{t}$ (33)
$\displaystyle\bm{B}_{i}+\bm{B}_{r}=\bm{B}_{t}+g\theta\bm{E}_{t}$ (34)

where the left-hand-side of the equals sign is the vacuum and the right-hand-
side is the DAI. In the vacuum $\bm{k}=\omega\hat{z}$ and in the DAI
$\bm{k}=k^{\pm}\hat{z}$ as above. The boundary condition for electric field
gives the simple condition $\bm{A}_{i}+\bm{A}_{r}=\bm{A}_{t}$. While the
boundary condition for magnetic field must be considered further. Substituting
back in for the vector potential

$\displaystyle\sum_{\alpha=\pm}-i\bm{k}_{i}\times\bm{A}_{i}^{\alpha}+i\bm{k}_{r}\times\bm{A}_{r}^{\alpha}=-i\bm{k}_{t}^{\alpha}\times\bm{A}_{t}^{\alpha}-ig\theta\omega\bm{A}_{t}^{\alpha}\quad\rightarrow$
(35)
$\displaystyle\sum_{\alpha}-i\omega\hat{z}\times\bm{A}_{i}^{\alpha}+i\omega\hat{z}\times\bm{A}_{r}^{\alpha}=-ik^{\alpha}\hat{z}\times\bm{A}_{t}^{\alpha}-ig\theta\omega\bm{A}_{t}^{\alpha}$
(36)

here let $\alpha=\pm$. Plugging in for the vector potentials in the chiral
basis and noting $\hat{z}\times\hat{e}_{\pm}=\pm i\hat{e}_{\pm}$

$\displaystyle\hat{e}_{+}:\omega I^{+}-\omega
R^{+}=(k^{+}-ig\theta\omega)T^{+}$ (37) $\displaystyle\hat{e}_{-}:\omega
I^{-}-\omega R^{-}=(k^{-}+ig\theta\omega)T^{-}$ (38)

Using $I^{+}+R^{+}=T^{+}$ and $I^{-}+R^{-}=T^{-}$ and solving for the $R^{+}$
gives:

$\displaystyle
R^{+}=\frac{\omega-(k^{+}-ig\theta\omega)}{\omega+(k^{+}-ig\theta\omega)}I^{+}$
(39)

Similarly for $R^{-}$:

$\displaystyle
R^{-}=\frac{-\omega+(k^{-}+ig\theta\omega)}{-\omega-(k^{-}+ig\theta\omega)}I^{-}$
(40)

recall $k^{\pm}=\omega\sqrt{\epsilon_{1}\pm\epsilon_{2}}$. The same can be
done for the transmission coefficients:

$\displaystyle T^{+}=\frac{2\omega}{\omega+(k^{+}-ig\theta\omega)}I^{+}$ (41)
$\displaystyle T^{-}=\frac{2\omega}{\omega+(k^{-}+ig\theta\omega)}I^{-}$ (42)

Finally, the Kerr rotation angle is related to the argument of the complex
reflection amplitudes$R^{\pm}=|R^{\pm}|e^{i\Delta^{\pm}}$ as

$\displaystyle\Theta_{K}=-\frac{1}{2}(\Delta^{+}-\Delta^{-})$ (43)

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Constraints imposed by $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. D 16, 1791 (1977b).
  * Chadha-Day  _et al._ [2022] F. Chadha-Day, J. Ellis, and D. J. E. Marsh, Axion dark matter: What is it and why now?, Science Advances 8, eabj3618 (2022), https://www.science.org/doi/pdf/10.1126/sciadv.abj3618 .
  * Nenno  _et al._ [2020] D. Nenno, C. Garcia, J. Gooth, C. Felser, and P. Narang, Axion physics in condensed-matter systems, Nat. Rev. Phys. 2, 682 (2020).
  * Li  _et al._ [2010] R. Li, J. Wang, X.-L. Qi, and S.-C. Zhang, Dynamical axion field in topological magnetic insulators, Nature Physics 6, 284 (2010), 0908.1537 .
  * Wang and Zhang [2013] Z. Wang and S.-C. Zhang, Chiral anomaly, charge density waves, and axion strings from Weyl semimetals, Physical Review B 87, 161107 (2013), 1207.5234 .
  * Roy and Sau [2015] B. Roy and J. D. Sau, Magnetic catalysis and axionic charge density wave in Weyl semimetals, Physical Review B 92, 125141 (2015), 1406.4501 .
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates, Physical Review B 83, 205101 (2011), 1007.0016 .
  * Pace  _et al._ [2021] S. D. Pace, C. Castelnovo, and C. R. Laumann, Dynamical Axions in $U(1)$ Quantum Spin Liquids, arXiv (2021), 2109.06890 .
  * Sekine and Nomura [2020] A. Sekine and K. Nomura, Axion electrodynamics in topological materials, Journal of Applied Physics 129, 141101 (2020).
  * Zyuzin  _et al._ [2012] A. A. Zyuzin, S. Wu, and A. A. Burkov, Weyl semimetal with broken time reversal and inversion symmetries, Physical Review B 85, 165110 (2012), 1201.3624 .
  * Essin  _et al._ [2009] A. M. Essin, J. E. Moore, and D. Vanderbilt, Magnetoelectric Polarizability and Axion Electrodynamics in Crystalline Insulators, Physical Review Letters 102, 146805 (2009), 0810.2998 .
  * Zyuzin and Burkov [2012] A. A. Zyuzin and A. A. Burkov, Topological response in weyl semimetals and the chiral anomaly, Phys. Rev. B 86, 115133 (2012).
  * Shi  _et al._ [2021] W. Shi, B. J. Wieder, H. L. Meyerheim, Y. Sun, Y. Zhang, Y. Li, L. Shen, Y. Qi, L. Yang, J. Jena, P. Werner, K. Koepernik, S. Parkin, Y. Chen, C. Felser, B. A. Bernevig, and Z. Wang, (A charge-density-wave topological semimetal, Nat. Phys. 17, 381 (2021).
  * Srivatsa and Ganesh [2018] N. S. Srivatsa and R. Ganesh, Excitonic collective modes in Weyl semimetals, Physical Review B 98, 165133 (2018), 1808.05233 .
  * Wieder  _et al._ [2020] B. J. Wieder, K.-S. Lin, and B. Bradlyn, Axionic band topology in inversion-symmetric Weyl-charge-density waves, Physical Review Research 2, 042010 (2020), 2004.11401 .
  * Sehayek  _et al._ [2020] D. Sehayek, M. Thakurathi, and A. A. Burkov, Charge density waves in Weyl semimetals, Physical Review B 102, 115159 (2020), 2007.07256 .
  * McKay and Bradlyn [2021] R. C. McKay and B. Bradlyn, Optical response from charge-density waves in Weyl semimetals, Physical Review B 104, 155120 (2021), 2106.09577 .
  * Curtis  _et al._ [2023] J. Curtis, I. Petrides, and P. Narang, Finite-Momentum Instability of Dynamical Axion Insulator, Phys. Rev. B 107, 205118 (2023).
  * Zhao  _et al._ [2020] B. Zhao, C. Guo, C. A. C. Garcia, P. Narang, and S. Fan, Axion-Field-Enabled Nonreciprocal Thermal Radiation in Weyl Semimetals, Nano Letters 20, 1923 (2020), 1911.08035 .
  * Otey  _et al._ [2010] C. R. Otey, W. T. Lau, S. Fan, _et al._ , Thermal rectification through vacuum, Physical review letters 104, 154301 (2010).
  * Gooth  _et al._ [2019] J. Gooth, B. Bradlyn, S. Honnali, C. Schindler, N. Kumar, J. Noky, Y. Qi, C. Shekhar, Y. Sun, Z. Wang, B. A. Bernevig, and C. Felser, Axionic charge-density wave in the weyl semimetal (TaSe4)2i, Nature 575, 315 (2019).
  * Yu  _et al._ [2021] J. Yu, B. J. Wieder, and C.-X. Liu, Dynamical piezomagnetic effect in time-reversal-invariant Weyl semimetals with axionic charge density waves, Physical Review B 104, 174406 (2021), 2008.10620 .
  * Zhu  _et al._ [2022] T. Zhu, H. Wang, H. Zhang, and D. Xing, Axionic surface wave in dynamical axion insulators, arXiv (2022), 2204.06332 .
  * Wu  _et al._ [2016] L. Wu, M. Salehi, N. Koirala, J. Moon, S. Oh, and N. P. Armitage, Quantized Faraday and Kerr rotation and axion electrodynamics of a 3D topological insulator, Science 354, 1124 (2016).
  * Ahn  _et al._ [2022] J. Ahn, S.-Y. Xu, and A. Vishwanath, Theory of optical axion electrodynamics and application to the Kerr effect in topological antiferromagnets, Nat. Comm. 13, 7614 (2022).
  * Yan and Felser [2017] B. Yan and C. Felser, Topological materials: Weyl semimetals, Annual Review of Condensed Matter Physics 8, 337 (2017), https://doi.org/10.1146/annurev-conmatphys-031016-025458 .
  * Burkov [2018] A. Burkov, Weyl metals, Annual Review of Condensed Matter Physics 9, 359 (2018), https://doi.org/10.1146/annurev-conmatphys-033117-054129 .
  * Note [1] In some cases it is also possible to have a separation of the Weyl nodes in energy of $Q_{0}$, leading to an additional contribution of $Q_{0}t$ in the static axion. This leads to additional topological responses but we will not consider this case here.
  * Zhang  _et al._ [2019] D. Zhang, M. Shi, T. Zhu, D. Xing, H. Zhang, and J. Wang, Topological axion states in the magnetic insulator ${\mathrm{mnbi}}_{2}{\mathrm{te}}_{4}$ with the quantized magnetoelectric effect, Phys. Rev. Lett. 122, 206401 (2019).
  * Kim  _et al._ [2021] S. Kim, R. McKay, _et al._ , Kramers-Weyl fermions in the chiral charge density wave material (TaSe4)2I (2021), arXiv:2108.10874 .
  * Nguyen  _et al._ [2022] Q. Nguyen, R. Duncan, _et al._ , Ultrafast x-ray scattering reveals composite amplitude collective mode in the Weyl charge density wave material (TaSe4)2I (2022), arXiv:2210.17483 .
  * Zhang  _et al._ [2022] Z. Zhang, F. Gao, _et al._ , Nonlinear coupled magnonics: Terahertz field-driven magnon upconversion (2022), arXiv:2207.07103 .
  * Zhang  _et al._ [2023] Z. Zhang, F. Gao, _et al._ , Three-wave mixing of anharmonically coupled magnons (2023), arXiv:2301.12555 .
  * Mootz  _et al._ [2022] M. Mootz, L. Luo, J. Wang, and I. Perakis, Visualization and quantum control of light-accelerated condensates by terahertz multi-dimensional coherent spectroscopy, Comm. Phys. 5, 47 (2022).
  * Rodriguez-Vega  _et al._ [2022] M. Rodriguez-Vega, Z.-X. Lin, A. Leonardo, A. Ernst, M. Vergniory, and G. Fiete, Light-Driven Topological and Magnetic Phase Transitions in Thin Layer Antiferromagnets, J. Phys. Chem. Lett. 14, 4152 (2022).
  * Peccei and Quinn [1977c] R. D. Peccei and H. R. Quinn, $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. Lett. 38, 1440 (1977c).
  * Schütte-Engel  _et al._ [2021] J. Schütte-Engel, D. J. Marsh, A. J. Millar, A. Sekine, F. Chadha-Day, S. Hoof, M. N. Ali, K. C. Fong, E. Hardy, and L. Š mejkal, Axion quasiparticles for axion dark matter detection, Journal of Cosmology and Astroparticle Physics 2021 (08), 066.
  * Chigusa  _et al._ [2021] S. Chigusa, T. Moroi, and K. Nakayama, Axion hidden-photon dark matter conversion into condensed matter axion, Journal of High Energy Physics 2021, 10.1007/jhep08(2021)074 (2021).