# Mendelian randomization with a binary exposure variable: interpretation and
presentation of causal estimates

Stephen Burgess 1,2 Corresponding author: Dr Stephen Burgess. Address: MRC
Biostatistics Unit, Cambridge Institute of Public Health, Robinson Way,
Cambridge, CB2 0SR, UK. Telephone: +44 1223 768259. Email:
<EMAIL_ADDRESS>Jeremy A Labrecque 3

1 MRC Biostatistics Unit, University of Cambridge, UK
2 Department of Public Health and Primary Care,
University of Cambridge, UK
3 Department of Epidemiology, Erasmus MC, Netherlands

Running head: Mendelian randomization with binary exposure

Keywords: Mendelian randomization, causal inference, instrumental variable,
effect estimation, power calculation.

Acknowledgements: The authors would like to thank Sonja A Swanson for her
contribution to an earlier draft of this manuscript.

## Abstract

Mendelian randomization uses genetic variants to make causal inferences about
a modifiable exposure. Subject to a genetic variant satisfying the
instrumental variable assumptions, an association between the variant and
outcome implies a causal effect of the exposure on the outcome. Complications
arise with a binary exposure that is a dichotomization of a continuous risk
factor (for example, hypertension is a dichotomization of blood pressure).
This can lead to violation of the exclusion restriction assumption: the
genetic variant can influence the outcome via the continuous risk factor even
if the binary exposure does not change. Provided the instrumental variable
assumptions are satisfied for the underlying continuous risk factor, causal
inferences for the binary exposure are valid for the continuous risk factor.
Causal estimates for the binary exposure assume the causal effect is a
stepwise function at the point of dichotomization. Even then, estimation
requires further parametric assumptions. Under monotonicity, the causal
estimate represents the average causal effect in ‘compliers’, individuals for
whom the binary exposure would be present if they have the genetic variant and
absent otherwise. Unlike in randomized trials, genetic compliers are unlikely
to be a large or representative subgroup of the population. Under homogeneity,
the causal effect of the exposure on the outcome is assumed constant in all
individuals; often an unrealistic assumption. We here provide methods for
causal estimation with a binary exposure (although subject to all the above
caveats). Mendelian randomization investigations with a dichotomized binary
exposure should be conceptualized in terms of an underlying continuous
variable.

Mendelian randomization is the use of genetic variants as instrumental
variables to test for or estimate the causal effect of a risk factor (referred
to here as an exposure) on an outcome using observational data (Davey Smith
and Ebrahim, 2003, Burgess and Thompson, 2015). The primary objective of
Mendelian randomization is to find modifiable exposures that are worthwhile
therapeutic targets and can be intervened on to improve health outcomes. An
instrumental variable must be associated with the exposure of interest
(relevance), only affects the outcome through the exposure (exclusion
restriction), and does not share any causes with the outcome
(exchangeability). Recently, several Mendelian randomization studies have
employed binary measures as the exposure variable. Examples include analyses
assessing the causal effect of cannabis initiation on schizophrenia (and of
schizophrenia on cannabis initiation) (Gage et al., 2017, Vaucher et al.,
2017), and of diabetes status on endometrial cancer (Nead et al., 2015). In
this short manuscript, we discuss issues relating to causal estimation in the
Mendelian randomization setting with a binary exposure. For ease of
presentation, we initially assume a single genetic variant is used as an
instrumental variable; this restriction is later relaxed.

The intended primary audience of this manuscript is Mendelian randomization
practitioners, and the aim of the manuscript is to communicate the practical
consequences of these methodological issues for Mendelian randomization
investigations. As such, we focus on methods and approaches that are likely to
be the most relevant to scenarios that are common in applied practice. In
particular, we focus on methods that can be performed using summarized data,
which comprise genetic associations with the exposure estimated using
regression methods, that are routinely reported by large consortia (Burgess et
al., 2013). Although our focus is on practitioners, we also provide technical
asides and references for methodologically-focused readers.

### Random assignment in a trial as a paradigm instrumental variable

Consider a double-blind, placebo-controlled randomized trial with two time-
fixed treatment arms (referred to as treatment and control) and complete
follow-up data. An intention-to-treat effect estimate is typically reported:
the causal effect of allocation to treatment as opposed to control. When there
is substantial non-compliance, investigators may be interested in testing
whether the treatment itself has an effect on the outcome (as opposed to
simply allocation to treatment), or in estimating the causal effect of the
treatment itself. Testing for a treatment or ‘per-protocol’ effect can be
achieved through the intention-to-treat analysis: unless random assignment
somehow affects the outcome directly (e.g., because blinding is broken or a
placebo effect is present), an association between treatment allocation and
the outcome will only arise if the treatment has a causal effect on the
outcome (Didelez and Sheehan, 2007). Estimating the average treatment effect
in the full study population further requires additional homogeneity
conditions (Hernán and Robins, 2006, Aronow and Carnegie, 2013, Wang and
Tchetgen Tchetgen, 2018); sufficient conditions are linearity of the
instrumental variable–exposure, instrumental variable–outcome and
exposure–outcome relationships with no effect heterogeneity. Without
additional conditions, only bounds for the average treatment effect are
obtainable (Balke and Pearl, 1997). These bounds can also be used to assess
the validity of a genetic variant as an instrumental variable (Ramsahai and
Lauritzen, 2011, Swanson et al., 2018), although this approach is rarely
informative in practice, and alternative ways of assessing instrument validity
(such as understanding the biological role of the genetic variant, and
assessing its associations with known confounders) are more likely to be
fruitful in practice (Burgess et al., 2015). Alternatively, investigators
often estimate an effect in a subgroup of the population under a weaker
assumption.

Specifically, we consider the subgroup of the population consisting of
‘compliers’ – individuals who would receive the treatment if allocated to
treatment, and would not receive treatment if allocated to not receive
treatment. The effect in this subgroup can be estimated under the assumption
that there are no defiers – individuals who would only take treatment if
randomly allocated not to do so, and who would not take treatment if allocated
to take it (Angrist et al., 1996). This is known as the monotonicity
assumption – allocation to taking the treatment can only increase the value of
the exposure, not decrease it. This effect, which can be estimated using
standard instrumental variable techniques, is known as the local average
treatment effect (LATE) or the complier average causal effect (CACE) in the
literature (Yau and Little, 2001). Of note, we cannot identify individual
compliers as we cannot see individuals’ treatment levels under both levels of
treatment allocation. However, it is possible to identify the proportion of
the study population who are compliers, and to describe relative
characteristics of the compliers compared to non-compliers using measured
baseline covariates (Angrist and Pischke, 2009). In well-designed randomized
trials, compliers are likely to be common, and the assumption that there are
no defiers is often considered reasonable.

### Who are the genetic ‘compliers’?

Monotonicity in the context of Mendelian randomization means that increasing
the number of variant alleles for an individual can only increase the exposure
from absent to present (or leave it constant), and can never decrease it. The
analogue of ‘compliers’ in Mendelian randomization are individuals who would
have the exposure present if they possess an exposure-increasing genetic
variant, but not otherwise. As genetic variants tend to have small effects on
phenotypic variables, such compliers are likely to be uncommon.

This means that the group of genetic compliers is not likely to be
representative of the general population. Also, the group of compliers may
well differ greatly between different study populations. As an example, folate
deficiency has been hypothesized as a causal risk factor for coronary heart
disease (Lewis et al., 2005). The complier population (and therefore the
instrumental variable estimate) would differ greatly in a population where
large numbers of people are borderline folate deficient compared with a
population where relatively few people are folate deficient. (A similar
problem would occur in randomized trials conducted in different populations.)
The analogous assumption in Mendelian randomization to the ‘no defiers’
assumption is that increases in the genotype variable would lead to increases
(or no change) in the exposure for all individuals in the population (or
equivalently, decreases or no change in the exposure for all individuals)
(Hernán and Robins, 2006). With a genetic variant that takes multiple values,
the equivalent assumption is that the exposure is a non-decreasing (or non-
increasing) function of the genetic variant. In this case (and in the case
with multiple genetic variants), the instrumental variable estimate is a
weighted average of LATEs (Angrist et al., 2000).

In the context of RCTs, even if individual compliers cannot be identified, the
subgroup of compliers may be of interest either because it represents a large
or representative subgroup of the population, or due to patterns of non-
compliance in the trial being anticipated to be repeated outside the trial
setting. However, in Mendelian randomization, the subgroup of genetic
‘compliers’ is unlikely to represent those individuals in the population who
would respond to a treatment that influences the target exposure, particularly
if the treatment has a greater effect on the risk factor than the genetic
variant. Hence, under the ‘no defiers’ assumption, the interpretation of a
causal estimate in a Mendelian randomization investigation in which the
instrumental variable assumptions are satisfied is that of an average causal
effect in those individuals whose exposure status would vary depending on
whether they have a particular genetic variant or not. We additionally note
that the subgroup of genetic compliers would differ between genetic variants.
This provides yet another reason why causal estimates based on different
genetic variants may vary even if all the genetic variants are valid
instruments.

### What is the true risk factor underlying the exposure?

The above interpretation assumes that the instrumental variable assumptions
are satisfied. These assumptions imply that the only influence of the
instrumental variable on the outcome is via the exposure – if the instrumental
variable changes, but the exposure stays the same, then the outcome should not
change. However, for most binary exposures used in Mendelian randomization
investigations, there is an underlying continuous risk factor for which the
binary variable is a dichotomization. As a simple example, the binary exposure
hypertension is a dichotomization of the continuous risk factor blood
pressure. In more complex examples, an underlying continuous latent variable
can be hypothesized even if it cannot be measured, such as a continuous
spectrum of sub-clinical mental health problems for the binary exposure
schizophrenia.

If the binary exposure is a dichotomization of a continuous risk factor, then
the instrumental variable assumptions are likely to be violated. For the
example of hypertension, if elevated blood pressure is a causal risk factor
for a particular outcome then genetic variants that are associated with blood
pressure will be associated with the outcome even in a population where no-one
suffers from clinically-defined hypertension. Hence, changes in the genetic
variants will lead to increases in blood pressure and consequently to changes
in the outcome even if the exposure status for hypertension remains fixed for
all individuals in the population. An instrumental variable for a continuous
exposure can only be an instrumental variable for the dichotomization of the
exposure if the exposure–outcome causal relationship is a strict stepwise
threshold at the point of dichotomization (in which case the dichotomized
exposure is a representation of the true risk factor). However, provided that
the instrumental variable assumptions are satisfied for the continuous risk
factor, testing for an association with the outcome is still a valid test of
the causal null hypothesis for the binary exposure.

There are two main consequences of this. First, such a Mendelian randomization
study should be conceptualized as an investigation into the (possibly latent)
underlying continuous risk factor, rather than the binary dichotomization of
this variable. At minimum, the instrumental variable assumptions should be
assessed with the continuous risk factor in mind. Second, a causal estimate
from a Mendelian randomization investigation with a dichotomized binary
exposure does not have a clear interpretation due to the binary exposure
variable not capturing the true causal relationship. There are several reasons
why a Mendelian randomization estimate may differ from the effect of an
intervention even for a continuous exposure (for example, genetic variants
have long-term influences acting from the beginning of life, whereas
interventions are more short-term and are applied to mature individuals)
(Burgess et al., 2012, Swanson et al., 2017). With a binary exposure, these
concerns are even greater.

### Causal estimation with a binary exposure

Despite this, suppose that we want to calculate a causal effect with a binary
exposure, under the assumption that the exposure has a stepwise effect on the
outcome. This may be because we truly believe in the homogeneity assumptions,
or we truly believe in the monotonicity assumption and regard the genetic
compliers as a worthwhile subgroup of the population in which to estimate an
average causal effect. Or, more likely, because a causal effect estimate is
required for pragmatic reasons, such as to perform a power calculation or to
inform policymakers of the expected impact of intervention on the exposure.
Other reasons for estimating a causal parameter include efficient testing of
the causal null hypothesis with multiple instrumental variables (under the
homogeneity assumptions, the two-stage least squares estimate, or equivalently
the inverse-variance weighted estimate, is the optimally efficient combination
of the instruments for testing for a causal effect (Wooldridge, 2009)) and use
of a robust method with multiple genetic variants (such as the MR-Egger method
(Bowden et al., 2015) or weighted median method (Bowden et al., 2016) – these
methods make weaker assumptions, not requiring all genetic variants to satisfy
the instrumental variable assumptions). If the binary exposure is a
dichotomization of a continuous risk factor, then power calculations are
likely to be conservative, as the effect of the genetic variant on the outcome
will not be fully captured by the binary exposure.

Two options for causal estimation are: i) estimating the effect on the outcome
per (say) 1% absolute increase in the probability of the exposure; ii)
estimating the effect on the outcome per (say) doubling of the probability (or
odds) of the exposure. We concentrate on estimation methods based on
regression (usually linear or logistic) for several reasons. First, often
researchers perform their analyses using summarized association estimates –
beta-coefficients from regression analyses of the exposure and outcome on a
genetic variant – and do not have access to individual-level data. These beta-
coefficients represent the average change in the trait (exposure or outcome)
per additional copy of the effect allele. Secondly, these approaches result in
causal estimates with a simple and relevant interpretation, and which can be
compared to estimates in the literature from other analytical approaches.
Thirdly, often there are technical restrictions on the data analysis – for
example, it may be necessary to fit a mixed model to account for relatedness
between individuals, to adjust for several principal components of ancestry,
or to provide a coordinated approach to analysis across different datasets.
These restrictions are easiest to accommodate in a regression framework. These
estimation procedures require strict linearity and homogeneity assumptions;
full details are available elsewhere (Hernán and Robins, 2006, Didelez and
Sheehan, 2007). The parametric assumptions for these two options are mutually
incompatible. Additionally, regression coefficients will generally be
variation dependent on the baseline risk, a nuisance parameter (Richardson et
al., 2017). If individual-level data are available, then alternative
approaches to estimation can be taken (Aronow and Carnegie, 2013, Wang and
Tchetgen Tchetgen, 2018).

If the genetic associations with the exposure are estimated using linear
regression, then they represent absolute changes in the prevalence of the
exposure. This enables estimation of the causal effect of an intervention in
the prevalence of the exposure on an absolute scale. It is sensible to scale
the causal effect to consider a modest increase in the prevalence of the
exposure (say a 1% or a 10% increase), as a unit increase would represent the
average causal effect of a population intervention from 0% prevalence of the
exposure to 100% prevalence – an unrealistic intervention in practice.
However, absolute associations with a binary variable do not make sense in
case-control settings (where cases are those with the exposure), as they
depend on the ratio of cases to controls chosen by the investigator.

If the genetic associations with the exposure are estimated using logistic
regression, then they represent log odds ratios. The causal estimate would
then represent the change in the outcome per unit change in the exposure on
the log odds scale. A unit increase in the log odds of a variable corresponds
to a 2.72 ($=\exp 1$)-fold multiplicative increase in the odds of the
variable. If the exposure is rare then the odds of the exposure is
approximately equal to the probability of the exposure. The causal estimate
represents the average change in the outcome per 2.72-fold increase in the
prevalence of the exposure (for example, an increase in the exposure
prevalence from 1% to 2.72%). It may be more interpretable to think instead
about the average change in the outcome per doubling (2-fold increase) in the
prevalence of the exposure. This can be obtained by multiplying the causal
estimate by 0.693 ($=\log_{e}2$).

### Discussion

In this short manuscript, we have discussed statistical issues for Mendelian
randomization with a binary exposure. A summary of the arguments made in the
paper is provided as Figure 1. Under the more plausible assumption of
monotonicity, the estimate from a Mendelian randomization study with a binary
exposure represents the average causal effect in ‘compliers’; the subgroup of
individuals for whom the presence or absence of the genetic variant used as an
instrument determines whether individuals have the exposure present or not.
Under the less plausible assumption of homogeneity, the estimate of the causal
effect only makes sense if the effect of the exposure on the outcome has a
strict stepwise form – only changes in whether the binary exposure is present
or absent will affect the outcome. If the binary exposure is a dichotomization
of a continuous variable, then the causal estimate does not have a clear
interpretation. In such a case, causal inferences will only be valid provided
that the instrumental variable assumptions are satisfied for the continuous
risk factor – in particular, if the effect of the genetic variant on the
outcome is completely mediated via the continuous risk factor. However, as the
effect of the genetic variant on the outcome is not completely mediated via
the binary exposure, power calculations are likely to be conservative.

In summary, applying Mendelian randomization with a binary exposure requires
careful consideration. When the binary exposure is a dichotomization of an
underlying continuous risk factor, causal assumptions should be assessed and
causal inferences should be conceptualized with respect to the underlying
continuous risk factor. Tests for causal effects may be achieved readily
without using the exposure information, but estimation procedures for a binary
exposure require strong assumptions that are unlikely to be biologically
plausible in common Mendelian randomization settings.

Funding: Stephen Burgess is supported by a Sir Henry Dale Fellowship jointly
funded by the Wellcome Trust and the Royal Society (Grant Number
204623/Z/16/Z).

Conflict of Interest: The authors declare that they have no conflict of
interest.

## References

  * Angrist et al. (2000) Angrist, J., K. Graddy, and G. Imbens (2000): “The interpretation of instrumental variables estimators in simultaneous equations models with an application to the demand for fish,” _Review of Economic Studies_ , 67, 499–527.
  * Angrist et al. (1996) Angrist, J., G. Imbens, and D. Rubin (1996): “Identification of causal effects using instrumental variables,” _Journal of the American Statistical Association_ , 91, 444–455.
  * Angrist and Pischke (2009) Angrist, J. and J. Pischke (2009): _Mostly harmless econometrics: an empiricist’s companion. Chapter 4: Instrumental variables in action: sometimes you get what you need_ , Princeton University Press.
  * Aronow and Carnegie (2013) Aronow, P. M. and A. Carnegie (2013): “Beyond LATE: Estimation of the average treatment effect with an instrumental variable,” _Political Analysis_ , 21, 492–506.
  * Balke and Pearl (1997) Balke, A. and J. Pearl (1997): “Bounds on treatment effects from studies with imperfect compliance,” _Journal of the American Statistical Association_ , 92, 1171–1176.
  * Bowden et al. (2015) Bowden, J., G. Davey Smith, and S. Burgess (2015): “Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression,” _International Journal of Epidemiology_ , 44, 512–525.
  * Bowden et al. (2016) Bowden, J., G. Davey Smith, P. C. Haycock, and S. Burgess (2016): “Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator,” _Genetic Epidemiology_ , 40, 304–314.
  * Brion et al. (2013) Brion, M.-J., K. Shakhbazov, and P. Visscher (2013): “Calculating statistical power in Mendelian randomization studies,” _International Journal of Epidemiology_ , 42, 1497–1501.
  * Burgess (2014) Burgess, S. (2014): “Sample size and power calculations in Mendelian randomization with a single instrumental variable and a binary outcome,” _International Journal of Epidemiology_ , 43, 922–929.
  * Burgess et al. (2012) Burgess, S., A. Butterworth, A. Malarstig, and S. Thompson (2012): “Use of Mendelian randomisation to assess potential benefit of clinical intervention,” _British Medical Journal_ , 345, e7325.
  * Burgess et al. (2013) Burgess, S., A. S. Butterworth, and S. G. Thompson (2013): “Mendelian randomization analysis with multiple genetic variants using summarized data,” _Genetic Epidemiology_ , 37, 658–665.
  * Burgess et al. (2015) Burgess, S., R. Scott, N. Timpson, G. Davey Smith, S. G. Thompson, and EPIC-InterAct Consortium (2015): “Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors,” _European Journal of Epidemiology_ , 30, 543–552.
  * Burgess and Thompson (2015) Burgess, S. and S. G. Thompson (2015): _Mendelian randomization: methods for using genetic variants in causal estimation_ , Chapman & Hall.
  * Davey Smith and Ebrahim (2003) Davey Smith, G. and S. Ebrahim (2003): “‘Mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?” _International Journal of Epidemiology_ , 32, 1–22.
  * Didelez and Sheehan (2007) Didelez, V. and N. Sheehan (2007): “Mendelian randomization as an instrumental variable approach to causal inference,” _Statistical Methods in Medical Research_ , 16, 309–330.
  * Gage et al. (2017) Gage, S. H., H. J. Jones, S. Burgess, J. Bowden, G. D. Smith, S. Zammit, and M. R. Munafò (2017): “Assessing causality in associations between cannabis use and schizophrenia risk: a two-sample Mendelian randomization study,” _Psychological Medicine_.
  * Hernán and Robins (2006) Hernán, M. and J. Robins (2006): “Instruments for causal inference: an epidemiologist’s dream?” _Epidemiology_ , 17, 360–372.
  * Lewis et al. (2005) Lewis, S., S. Ebrahim, and G. Davey Smith (2005): “Meta-analysis of MTHFR 677C - T polymorphism and coronary heart disease: does totality of evidence support causal role for homocysteine and preventive potential of folate?” _British Medical Journal_ , 331, 1053.
  * Nead et al. (2015) Nead, K. T., S. J. Sharp, D. J. Thompson, J. N. Painter, D. B. Savage, R. K. Semple, A. Barker, J. R. Perry, J. Attia, A. M. Dunning, et al. (2015): “Evidence of a causal association between insulinemia and endometrial cancer: a Mendelian randomization analysis,” _Journal of the National Cancer Institute_ , 107, djv178.
  * Ramsahai and Lauritzen (2011) Ramsahai, R. and S. Lauritzen (2011): “Likelihood analysis of the binary instrumental variable model,” _Biometrika_ , 98, 987–994.
  * Richardson et al. (2017) Richardson, T. S., J. M. Robins, and L. Wang (2017): “On modeling and estimation for the relative risk and risk difference,” _Journal of the American Statistical Association_ , 112, 1121–1130.
  * Swanson et al. (2018) Swanson, S. et al. (2018): “Partial identification of the average treatment effect using instrumental variables: review of methods for binary instruments, treatments, and outcomes,” _Journal of the American Statistical Association_.
  * Swanson et al. (2017) Swanson, S. A., H. Tiemeier, M. A. Ikram, and M. A. Hernán (2017): “Nature as a trialist?: Deconstructing the analogy between Mendelian randomization and randomized trials,” _Epidemiology_ , 28, 653–659.
  * Vaucher et al. (2017) Vaucher, J., B. J. Keating, A. M. Lasserre, W. Gan, D. Lyall, J. Ward, D. J. Smith, J. Pell, N. Sattar, G. Pare, and M. Holmes (2017): “Cannabis use and risk of schizophrenia: a mendelian randomization study,” _Molecular Psychiatry_.
  * Wang and Tchetgen Tchetgen (2018) Wang, L. and E. Tchetgen Tchetgen (2018): “Bounded, efficient and triply robust estimation of average treatment effects using instrumental variables,” _arXiv_ , 1611.09925.
  * Wooldridge (2009) Wooldridge, J. (2009): _Introductory econometrics: A modern approach. Chapter 15: Instrumental variables estimation and two stage least squares_ , South-Western, Nashville, TN.
  * Yau and Little (2001) Yau, L. H. and R. J. Little (2001): “Inference for the complier-average causal effect from longitudinal data subject to noncompliance and missing data, with application to a job training assessment for the unemployed,” _Journal of the American Statistical Association_ , 96, 1232–1244.

Figure 1: Flow diagram illustrating the steps needed to consider when
considering whether to estimate a parameter in a Mendelian randomization
investigation or not with a binary risk factor.


#  A Liouville type theorem for ancient Lagrangian Mean Curvature flows

Arunima Bhattacharya, Micah Warren, and Daniel Weser  Department of
Mathematics, Phillips Hall
the University of North Carolina at Chapel Hill, NC<EMAIL_ADDRESS>Department of Mathematics
University of Oregon, Eugene, OR 97403<EMAIL_ADDRESS>Department of
Mathematics, Phillips Hall
the University of North Carolina at Chapel Hill, NC<EMAIL_ADDRESS>
###### Abstract.

We prove a Liouville type result for convex solutions of the Lagrangian mean
curvature flow with restricted quadratic growth assumptions at antiquity on
the solutions.

##  1\. Introduction

A family of Lagrangian submanifolds
$X(x,t):\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{C}^{n}$ is said to evolve by
mean curvature if

$(X_{t})^{\bot}=\Delta_{g}X=\vec{H},$

where $\vec{H}$ denotes the mean curvature vector of the Lagrangian
submanifold. After a change of coordinates, one can locally write
$X(x,t)=(x,Du(x,t))$ such that $\Delta_{g}X=J\nabla_{g}\Theta$ (see [HL82]):
Here $g=I_{n}+(D^{2}u)^{2}$ is the induced metric on $(x,Du(x))$, $J$ is the
almost complex structure on $\mathbb{C}^{n}$, and $\Theta$ is the Lagrangian
angle given by

$\Theta=\sum_{i=1}^{n}\arctan\lambda_{i},$

where $\lambda_{i}$ are the eigenvalues of the Hessian $D^{2}u$. This results
in a local potential $u(x,t)$ evolving by the parabolic equation

(1.1) $\displaystyle u_{t}$ $\displaystyle=\sum_{i=1}^{n}\arctan\lambda_{i}.$

We are free to add any function that does not depend on $x$ to the right-hand
side, since this does not change the flow of the gradient graph of $u$. Thus,
we will consider the equation

(1.2) $u_{t}=\sum_{i=1}^{n}\arctan\lambda_{i}-\Theta_{0}$

for a fixed $\Theta_{0}$. In particular, under this flow, solutions to the
special Lagrangian equation

(1.3) $\sum_{i=1}^{n}\arctan\lambda_{i}=\Theta_{0}$

will be stationary.

Our main result in this paper is the following:

######  Theorem 1.1.

Let $u$ be an entire ancient solution to (1.2). Suppose that $u$ is convex and
satisfies the following growth condition at antiquity: there exists an $R_{0}$
such that

$\lim\sup_{t\rightarrow-\infty}\sup_{x\in\mathbb{R}^{n}}\frac{\left|u(x,t)\right|}{\left|x\right|^{2}+R_{0}}<\frac{1}{(6\sqrt{n}+2)^{2}}.$

Then $u$ is a quadratic polynomial

In essence, this theorem states that if a convex ancient solution has the
growth of a small quadratic polynomial for times far enough in the past, then
the solution must in fact be a quadratic polynomial. Although the growth
condition is somewhat restrictive, it applies to the scalar potential $u$, not
to the gradient graph, and requires only convexity without additional
constraints on the Hessian.

We derive a Hessian estimate for solutions of (1.1) via a pointwise approach
assuming a quadratic growth condition at antiquity. The method used here
follows a parabolic adaptation of Korevaar’s maximum principle method as in
Evan-Spruck [ES92]. A similar maximum principle estimate for the elliptic
special Lagrangian equation was carried out in Warren-Yuan [WY08], inspired by
Korevaar’s [Kor86] estimate for minimal surface equations. As in [WY08], our
method requires a smallness condition on the growth. Once a Hessian bound is
established, Hölder regularity and the Liouville property follow from Nguyen-
Yuan [NY11].

In 1968, Bombieri, De Giorgi, and Miranda [BDGM69] established an interior
estimate bound for classical solutions of the minimal surface equation in
$\mathbb{R}^{n}$. In the 1970s, Trudinger [Tru72] provided a new and simpler
derivation of this estimate and partly developed in the process some new
techniques applicable to the study of hypersurfaces in general. In the 1980s,
Korevaar [Kor86] found a strikingly simple pointwise argument. Korevaar’s
technique was adapted to prove estimates for the parabolic version, namely
mean curvature flow, by Ecker-Huisken [EH91], Evans-Spruck [ES92], Colding-
Minicozzi [CMI04] in codimension one.

Mean curvature flow in higher codimensions including the Lagrangian mean
curvature flow has been studied by many authors [CCY12, CCY13, SW02, Wan02,
TTW23, Wan24, LSS24]. However, in higher codimensions, an adaption of
Korevaar’s method with no restriction on the height remains elusive, as the
diagonization step is difficult to overcome. To diagonalize the differential
of a vector-valued function by rotating the base, one requires that the
gradients of each of the functions in the vector be mutually orthogonal, which
fails in general. Interior Hessian estimates for the elliptic special
Lagrangian equation and the variable phase Lagrangian mean curvature equation
have been established via an integral approach in [BC03, BCGJ03, WY09a, WY09b,
CWY09, WY10, WY14, CSY19, Bha21, Bha22, BMS24, Zho22, BW24], via a compactness
approach in [Li19], via a doubling approach in [Sha24]. The Liouville result
of [WY08] is generalized to remove the small quadratic constraint via a
compactness method in [Din23]. A pointwise approach to proving interior
Hessian estimates without any constraints on the growth remains an open
question even for the elliptic special Lagrangian equation.

Acknowledgments. AB acknowledges the support of the Simons Foundation grant
MP-TSM-00002933. DW acknowledges the support of the NSF RTG DMS-2135998 grant.

##  2\. Jacobi Inequality

Let $g$ denote the induced metric of the embedding $X$, which in the graphical
case takes the form $g=I_{n}+D^{2}u(\cdot,t)D^{2}u(\cdot,t)$. If we take a
derivative of (1.2), we find

$u_{tk}=g^{ij}u_{ijk},$

so that, if we define the operator

$L:=\partial_{t}-g^{ij}\partial_{ij},$

we see that equivalently

(2.1) $Lu_{k}=0.$

We will let $V=\sqrt{|\det g|}$ denote the volume element, and we will define
$b=V^{1/n}$ to be the volume element raised to the power of $1/n$. In the
following lemma, we show that a Jacobi inequality holds along the “heat flow”
of $b$ with respect to $L$, which will be a crucial ingredient in the proof of
the main theorem.

######  Lemma 2.1.

For a convex solution of (1.1), the volume element $b=(\sqrt{|\det g|})^{1/n}$
satisfies

(2.2) $Lb+2\frac{|\nabla_{g}b|^{2}}{b}\leq 0.$

###### Proof.

We directly compute

(2.3) $\displaystyle\partial_{t}b$
$\displaystyle=\frac{1}{2n}\mathrm{tr}(g^{-1}\partial_{t}g)b$ (2.4)
$\displaystyle\partial_{i}b$
$\displaystyle=\frac{1}{2n}\mathrm{tr}(g^{-1}\partial_{i}g)b$ (2.5)
$\displaystyle\partial_{ij}b$
$\displaystyle=\frac{1}{2n}\mathrm{tr}(g^{-1}\partial_{ij}g)b+\frac{1}{2n}\mathrm{tr}(\partial_{j}g^{-1}\partial_{i}g)b+\frac{\partial_{i}b\partial_{j}b}{b}.$

Therefore, we find the initial expansion

(2.6)
$Lb=\frac{1}{2n}\mathrm{tr}(g^{-1}Lg)b-\frac{1}{2n}g^{ij}\mathrm{tr}(\partial_{j}g^{-1}\partial_{i}g)b-\frac{|\nabla_{g}b|^{2}}{b}.$

Now we will compute $Lg$. To do so, we fix a point $p$ and rotate coordinates
so that $D^{2}u$ is diagonalized at $p$. (Note; a rotation in the domain of a
scalar function generates a transpose rotation of the gradient of that
function. The conjugation by these is what diagonalizes the Hessian.) We
denote by $\lambda_{k}$ the eigenvalues of $D^{2}u(p)$, so that
$D^{2}u(p)=\textrm{diag}(\lambda_{1},\ldots,\lambda_{n})$. Note that
$\lambda_{k}\geq 0$ by our convexity hypothesis. In these coordinates at the
point $p$, the metric $g$ takes the form

(2.7) $g_{kl}=\delta_{kl}+Du_{k}\cdot Du_{l}=(1+\lambda_{k}^{2})\delta_{kl},$

and thus

(2.8) $g^{kl}=(1+\lambda_{k}^{2})^{-1}\delta^{kl}.$

Applying these identities at the point $p$, one obtains

(2.9) $\displaystyle\partial_{t}g_{kl}$
$\displaystyle=(\lambda_{k}+\lambda_{l})\partial_{t}u_{kl}$ (2.10)
$\displaystyle\partial_{ij}g_{kl}$
$\displaystyle=(\lambda_{k}+\lambda_{l})u_{ijkl}+Du_{ik}\cdot
Du_{jl}+Du_{jk}\cdot Du_{il},$

so that

(2.11) $\displaystyle Lg_{kl}$
$\displaystyle=(\lambda_{k}+\lambda_{l})\partial_{t}u_{kl}-g^{ij}\Big{(}(\lambda_{k}+\lambda_{l})u_{ijkl}+Du_{ik}\cdot
Du_{jl}+Du_{jk}\cdot Du_{il}\Big{)}$ (2.12)
$\displaystyle=(\lambda_{k}+\lambda_{l})\partial_{t}u_{kl}-g^{ii}(\lambda_{k}+\lambda_{l})u_{iikl}+2g^{ii}Du_{ik}\cdot
Du_{il}.$

Thus, at the point $p$ we find

(2.13) $\displaystyle\mathrm{tr}(g^{-1}Lg)$ $\displaystyle=g^{kk}Lg_{kk}$
(2.14)
$\displaystyle=2g^{kk}\lambda_{k}\partial_{t}u_{kk}-2g^{ii}g^{kk}\lambda_{k}u_{iikk}-2g^{ii}g^{kk}|Du_{ik}|^{2}$
(2.15) $\displaystyle=2g^{kk}\lambda_{k}Lu_{kk}-2g^{ii}g^{kk}|Du_{ik}|^{2}.$

Now we compute $Lu_{kk}$. Differentiating $Lu_{k}=0$, we obtain

(2.16) $Lu_{kk}=(\partial_{k}g^{ij})u_{ijk}.$

Differentiating the identity $\delta_{a}^{c}=g_{ab}g^{bc}$, one obtains at the
point $p$ that

(2.17) $\partial_{k}g^{ij}=-g^{ii}g^{jj}(\lambda_{i}+\lambda_{j})u_{ijk},$

so that

(2.18) $\displaystyle Lu_{kk}$ $\displaystyle=(\partial_{k}g^{ij})u_{ijk}$
(2.19) $\displaystyle=-g^{ii}g^{jj}(\lambda_{i}+\lambda_{j})u_{ijk}^{2}.$

Combining (2.15) and (2.19), we find at the point $p$ that

(2.20) $\displaystyle\mathrm{tr}(g^{-1}Lg)$
$\displaystyle=2g^{kk}\lambda_{k}Lu_{kk}-2g^{ii}g^{kk}|Du_{ik}|^{2}$ (2.21)
$\displaystyle=-2g^{ii}g^{jj}g^{kk}(\lambda_{i}+\lambda_{j})\lambda_{k}Lu_{ijk}^{2}-2g^{ii}g^{kk}|Du_{ik}|^{2}.$

Next, using (2.17), we directly compute

(2.22) $\displaystyle g^{ij}\mathrm{tr}(\partial_{j}g^{-1}\partial_{i}g)$
$\displaystyle=g^{ii}\partial_{i}g^{kl}\partial_{i}g_{kl}$ (2.23)
$\displaystyle=-g^{ii}g^{kk}g^{ll}(\lambda_{k}+\lambda_{l})^{2}u_{ikl}^{2},$

so that, inserting (2.21) and (2.23) into (2.6), we obtain

(2.24) $\displaystyle Lb$
$\displaystyle=\frac{1}{2n}\mathrm{tr}(g^{-1}Lg)b-\frac{1}{2n}g^{ij}\mathrm{tr}(\partial_{j}g^{-1}\partial_{i}g)b-\frac{|\nabla_{g}b|^{2}}{b}$
(2.25)
$\displaystyle=\frac{b}{2n}\bigg{(}-2g^{ii}g^{jj}g^{kk}(\lambda_{i}+\lambda_{j})\lambda_{k}u_{ijk}^{2}-2g^{ii}g^{kk}|Du_{ik}|^{2}$
(2.26) $\displaystyle\hskip
56.9055pt+g^{ii}g^{kk}g^{ll}(\lambda_{k}+\lambda_{l})^{2}u_{ikl}^{2}\bigg{)}-\frac{|\nabla_{g}b|^{2}}{b}$
(2.27)
$\displaystyle=\frac{b}{2n}\bigg{(}-2g^{ii}g^{jj}g^{kk}(\lambda_{i}+\lambda_{j})\lambda_{k}u_{ijk}^{2}-2g^{ii}\delta^{jj}g^{kk}u_{ijk}^{2}$
(2.28) $\displaystyle\hskip
56.9055pt+g^{ii}g^{kk}g^{ll}(\lambda_{k}+\lambda_{l})^{2}u_{ikl}^{2}\bigg{)}-\frac{|\nabla_{g}b|^{2}}{b}$
(2.29)
$\displaystyle=\frac{b}{2n}\bigg{(}-2g^{ii}g^{jj}g^{kk}(\lambda_{i}+\lambda_{j})\lambda_{k}u_{ijk}^{2}-2g^{ii}g^{jj}g^{kk}(1+\lambda_{j}^{2})u_{ijk}^{2}$
(2.30) $\displaystyle\hskip
56.9055pt+g^{ii}g^{kk}g^{ll}(\lambda_{k}+\lambda_{l})^{2}u_{ikl}^{2}\bigg{)}-\frac{|\nabla_{g}b|^{2}}{b}$
(2.31)
$\displaystyle=\frac{b}{2n}\bigg{(}-4g^{ii}g^{jj}g^{kk}\lambda_{i}\lambda_{k}u_{ijk}^{2}-2g^{ii}g^{jj}g^{kk}\lambda_{j}^{2}u_{ijk}^{2}-2g^{ii}g^{jj}g^{kk}u_{ijk}^{2}$
(2.32) $\displaystyle\hskip
56.9055pt+2g^{ii}g^{kk}g^{ll}\lambda_{k}^{2}u_{ikl}^{2}+2g^{ii}g^{kk}g^{ll}\lambda_{k}\lambda_{l}u_{ikl}^{2}\bigg{)}-\frac{|\nabla_{g}b|^{2}}{b}$
(2.33)
$\displaystyle=\frac{b}{2n}\bigg{(}-2g^{ii}g^{jj}g^{kk}\lambda_{i}\lambda_{k}u_{ijk}^{2}-2g^{ii}g^{jj}g^{kk}u_{ijk}^{2}\bigg{)}-\frac{|\nabla_{g}b|^{2}}{b}.$

Finally, we use Cauchy-Schwarz to see that

(2.34)
$\displaystyle\frac{|\nabla_{g}b|^{2}}{b}=\frac{b}{n^{2}}g^{ii}\Big{(}g^{kk}\lambda_{k}u_{ikk}\Big{)}^{2}\leq\frac{b}{n}g^{ii}g^{kk}g^{kk}\lambda_{k}^{2}u_{ikk}^{2},$

which, by comparing with (2.33) and using the fact that $\lambda_{k}\geq 0$ by
convexity, concludes the poof of the lemma. ∎

##  3\. Oscillation and gradient bounds forward in time

In this section, we will prove interior oscillation and gradient bounds
forward in time for convex solutions of (1.2) defined on parabolic cylinders
of the form

$\overline{B}_{R}(0)\times[0,\frac{1}{n}]\subset\mathbb{R}^{n}\times\mathbb{R}.$

First, we prove an interior height bound forward in time for general solutions
of (1.1).

######  Proposition 3.1.

Suppose that $u$ is a solution of (1.1) on
$\overline{B}_{R}(0)\times[0,\frac{1}{n}]$. Then we have

(3.1)
$u(0,\frac{1}{n})\leq\arctan\left(\frac{\pi}{R^{2}}\right)+\max_{\overline{B}_{R}(0)\times\left\\{0\right\\}}u\left(x,0\right).$

###### Proof.

Up to subtracting a constant, we will assume that $\max u=0$ at $t=0.$ Notice
that by simply integrating equation (1.1) at each point on $\partial B_{R}(0)$
we may obtain the bound

(3.2) $u(x,t)\leq\frac{tn\pi}{2}\text{ on }\partial
B_{R}(0)\times[0,\frac{1}{n}].$

We will construct a supersolution to obtain the interior bound. Let

$w=\frac{tn\pi}{2}\left(\frac{\left|x\right|}{R}\right)^{2}+nt\arctan\left(\frac{tn\pi}{R^{2}}\right),$

which satisfies

$w_{t}-\sum\arctan\lambda_{i}\geq 0.$

By (3.2), $w\geq u$ on the boundary, and hence also in the interior of
$\overline{B}_{R}(0)\times[0,\frac{1}{n}]$. Thus

$\displaystyle u(0,\frac{1}{n})$ $\displaystyle\leq
w(0,\frac{1}{n})\leq\arctan\left(\frac{\pi}{R^{2}}\right).$

∎

######  Corollary 3.1.

Suppose for any $\Theta_{0}$ the function $u$ is a convex solution of (1.2) on

$\overline{B}_{2R+1}(0)\times[0,\frac{1}{n}]\subset\mathbb{R}^{n}\times\mathbb{R}.$

Then,

$\max_{\overline{B}_{1}(0)\times[0,\frac{1}{n}]}|Du|\leq\frac{1}{R}\left[M+\arctan\left(\frac{\pi}{R^{2}}\right)\right]$

where $M$ is the oscillation of $u$ at $t=0.$

###### Proof.

By subtracting a constant, we will assume that
$\max_{\overline{B}_{2R+1}(0)\times[0,T)}u=0$, so that
$\min_{\overline{B}_{2R+1}(0)\times[0,T)}u=-M.$ Adding $\Theta_{0}t$, we may
also assume that $u$ is a solution of (1.1). From (3.1), we have

$u(x,t)\leq\arctan\left(\frac{\pi}{R^{2}}\right)\text{ on
}B_{R+1}(0)\times[0,\frac{1}{n}].$

Since $u$ is convex, we know $\Theta\geq 0$, and hence $u$ can only increase
along the flow of type (1.1). Therefore, we see that the oscillation is no
more than $M+\arctan\left(\frac{\pi}{R^{2}}\right)$ on
$B_{R+1}(0)\times[0,\frac{1}{n}].$ The conclusion follows by applying the
convexity condition to any point on $B_{1}(0)$, noting that the derivative
bound from convexity depends only on the oscillation of each time slice.
However, since solutions of (1.1) and (1.2) only differ by a constant in each
time slice, we also obtain the conclusion for solutions of (1.2). ∎

##  4\. Hessian bound

In this section, we prove an interior Hessian bound forward in time for convex
solutions of (1.1) and (1.2). The first proposition applies to convex
solutions with sufficiently small gradient $|Du|$, while the second applies to
convex solutions with bounded oscillation.

######  Proposition 4.1.

Let $u$ be a convex solution of (1.1) or (1.2) on
$B_{1}(0)\times[0,\frac{1}{n}]$. Suppose that for some $\alpha,\gamma>0$ we
have

$\alpha>\frac{3n}{2},$ $\alpha\gamma<1,$

and

$\left|Du\right|^{2}<\gamma\text{ on }B_{1}(0)\times[0,\frac{1}{n}].$

Then for any $K>0$

$\lambda_{\max}^{2}(0,\frac{1}{n})\leq\frac{e^{2nK}\left(\frac{2\alpha}{2\alpha-3n}\right)^{n}}{\left(e^{K\left(1-\alpha\gamma\right)}-1\right)^{2n}}.$

###### Proof.

Let

$\displaystyle b$ $\displaystyle=V^{1/n}$ $\displaystyle\phi$
$\displaystyle=\left[\alpha\left|Du\right|^{2}-\alpha\gamma+nt(1-\left|x\right|^{2})\right]^{+}$
$\displaystyle f$ $\displaystyle=e^{K\phi}-1$ $\displaystyle\eta$
$\displaystyle=f\circ\phi,$

and consider the function

$h=\eta b.$

Notice that

$-\alpha\gamma\leq\alpha\left|Du\right|^{2}-\alpha\gamma\leq 0,$

so $\phi$ is 0 when $t=0$ (and thus so is $\eta$). Additionally, $\phi$
vanishes on the boundary $\partial B_{1}(0)$ for any time.

Now, since $-\alpha\gamma+1>0$, the function $\eta$ will be positive
somewhere, and we let $\left(x_{0},t_{0}\right)$ be a point where a positive
maximum for $h$ occurs on $B_{1}(0)\times[0,\frac{1}{n}]$. Then we must have

$Lh=h_{t}-g^{ij}h_{ij}\geq 0.$

From spatial maximality, we obtain

$\eta_{i}=-\frac{\eta b_{i}}{b},$

which implies

$\displaystyle Lh$ $\displaystyle=\eta_{t}b+\eta
b_{t}-g^{ij}\eta_{ij}b-2g^{ij}\eta_{i}b_{j}-\eta g^{ij}b_{ij}$
$\displaystyle=\left(\eta_{t}-g^{ij}\eta_{ij}\right)b+\eta\left(b_{t}+2g^{ij}\frac{b_{i}}{b}b_{j}-g^{ij}b_{ij}\right).$

By (2.2), we know that

$\left(b_{t}+2g^{ij}\frac{b_{i}}{b}b_{j}-g^{ij}b_{ij}\right)\leq 0,$

so we conclude at $\left(x_{0},t_{0}\right)$ that

$\eta_{t}-g^{ij}\eta_{ij}\geq 0.$

We compute

$\displaystyle\eta_{i}$
$\displaystyle=Ke^{K\phi}\left(2\alpha\sum_{k=1}^{n}u_{k}u_{ki}-2ntx_{i}\right)$
$\displaystyle\eta_{ij}$
$\displaystyle=Ke^{K\phi}\left(2\alpha\sum_{k=1}^{n}u_{kj}u_{ki}+2\alpha\sum_{k=1}^{n}u_{kij}u_{k}-2nt\delta_{ij}\right)$
$\displaystyle\hskip
14.22636pt+K^{2}e^{K\phi}\left(2\alpha\sum_{k=1}^{n}u_{k}u_{ki}-2ntx_{i}\right)\left(2\alpha\sum_{k=1}^{n}u_{k}u_{kj}-2ntx_{j}\right)$
$\displaystyle\eta_{t}$
$\displaystyle=Ke^{K\phi}\left(2\alpha\sum_{k=1}^{n}u_{k}u_{kt}+n\left(1-\left|x\right|^{2}\right)\right),$

so at $\left(x_{0},t_{0}\right)$ we have

$\displaystyle 0\leq\frac{\eta_{t}-g^{ij}\eta_{ij}}{Ke^{K\phi}}$
$\displaystyle=\left(2\alpha\sum_{k=1}^{n}u_{k}u_{kt}+n\left(1-\left|x\right|^{2}\right)\right)$
$\displaystyle\hskip
21.33955pt-g^{ij}\left(2\alpha\sum_{k=1}^{n}u_{kj}u_{ki}+2\alpha\sum_{k=1}^{n}u_{kij}u_{k}-2nt\delta_{ij}\right)$
$\displaystyle\hskip
21.33955pt-g^{ij}K\left(2\alpha\sum_{k=1}^{n}u_{k}u_{ki}-2ntx_{i}\right)\left(2\alpha\sum_{k=1}^{n}u_{k}u_{kj}-2ntx_{j}\right).$

Using (2.1), discarding the final term, and diagonalizing the Hessian (which
will change the direction of the gradient at $(x_{0},t_{0})$ but not its
magnitude) we obtain

$\displaystyle 0$ $\displaystyle\leq
n\left(1-\left|x\right|^{2}\right)-2\alpha\sum\frac{\lambda_{i}^{2}}{1+\lambda_{i}^{2}}+2nt\sum\frac{1}{1+\lambda_{i}^{2}},$

which can be written as

$2\alpha\sum\frac{\lambda_{i}^{2}}{1+\lambda_{i}^{2}}\leq\ 3n.$

Therefore, we conclude at $\left(x_{0},t_{0}\right)$ that any eigenvalue
$\lambda_{i}$ satisfies the bound

$\frac{\lambda_{i}^{2}}{1+\lambda_{i}^{2}}\leq\frac{3n}{2\alpha},$

which is equivalent to

$\lambda_{i}^{2}\leq\frac{3n}{2\alpha-3n}.$

Thus,

$\displaystyle V\left(x_{0},t_{0}\right)$
$\displaystyle=\prod\left(1+\lambda_{i}^{2}\right)^{1/2}$
$\displaystyle\leq\left[1+\left(\frac{3n}{2\alpha-3n}\right)\right]^{n/2}$
$\displaystyle=\left(\frac{2\alpha}{2\alpha-3n}\right)^{n/2},$

and so

$b\left(x_{0},t_{0}\right)\leq\left(\frac{2\alpha}{2\alpha-3n}\right)^{1/2}.$

Now, for all $\left(x,t\right)$ we have $\phi\leq 1$ and hence

$h\left(x_{0},t_{0}\right)\leq
e^{K}\left(\frac{2\alpha}{2\alpha-3n}\right)^{1/2},$

from which it follows that

$h(0,T)\leq e^{K}\left(\frac{2\alpha}{2\alpha-3n}\right)^{1/2}.$

Finally, observe that

$\phi(0,\frac{1}{n})\geq-\alpha\gamma+1,$

which yields

$b(0,\frac{1}{n})\leq\frac{e^{K}\left(\frac{2\alpha}{2\alpha-3n}\right)^{1/2}}{e^{K\left(1-\alpha\gamma\right)}-1}.$

Therefore, the maximum eigenvalue $\lambda_{\max}$ satisfies

$1+\lambda_{\max}^{2}\leq\frac{e^{2nK}\left(\frac{2\alpha}{2\alpha-3n}\right)^{n}}{\left(e^{K\left(1-\alpha\gamma\right)}-1\right)^{2n}}.$

∎

We now combine Proposition 4.1 with Corollary 3.1 to obtain a Hessian bound
from the oscillation.

######  Proposition 4.2.

Suppose that $u$ is a convex solution of (1.2) on
$B_{2R+1}(0)\times[0,\frac{1}{n}]$ with

$R=3\sqrt{n}$

and $u$ satisfies

$-1\leq u\leq 1$

on $B_{2R+1}\times\left\\{0\right\\}$. Then

$D^{2}u(0,\frac{1}{n})\leq C(n).$

###### Proof.

By Corollary 3.1, we have the bound

$\max_{B_{1}(0)\times[0,\frac{1}{n}]}|Du|\leq\frac{1}{3\sqrt{n}}\left[2+\arctan\left(\frac{\pi}{9n}\right)\right].$

We now apply the Hessian bound of Proposition 4.1 with

$\displaystyle T$ $\displaystyle=\frac{1}{n}$ $\displaystyle\gamma$
$\displaystyle=\frac{1}{9n}\left(2+\arctan\left(\frac{\pi}{9n}\right)\right)^{2}<\frac{0.61}{n}$
$\displaystyle\alpha$ $\displaystyle=1.6n$ $\displaystyle K$ $\displaystyle=1$

to conclude

$\lambda_{\max}^{2}(0,\frac{1}{n})\leq\frac{e\left(16\right)^{n}}{\left(e^{\left(1-\alpha\gamma\right)}-1\right)^{2n}}=C(n).$

∎

##  5\. Proof of the Liouville Theorem

To prove our main theorem we will show that for any $\left(x_{0},t_{0}\right)$
we have $D^{2}u(x_{0},t_{0})\leq C(n).$ Without loss of generality we may
choose $t_{0}=0.$ For fixed $x_{0}$, consider

$\tilde{u}_{\lambda}(x,t)=\frac{1}{\lambda^{2}}u\left(\lambda(x-x_{0}),\lambda^{2}t\right),$

which restricts to a solution of (1.1) on
$B_{6\sqrt{n}+1}(0)\times[-\frac{1}{n},0]$ with

(5.1) $\sup_{x\in
B_{6\sqrt{n}+1}(0)}\left|\tilde{u}_{\lambda}(x,-\frac{1}{n})\right|\leq\frac{1}{\lambda^{2}}\sup_{z\in
B_{\left(6\sqrt{n}+1\right)\left(\lambda+\left|x_{0}\right|\right)}(0)}\left|u\left(z,-\lambda^{2}\frac{1}{n}\right)\right|.$

By the growth at antiquity condition, for $\lambda$ large enough and
$\varepsilon_{0}<1$ small enough we have

$\displaystyle\sup_{z\in
B_{\left(6\sqrt{n}+1\right)\left(\lambda+\left|x_{0}\right|\right)}(0)}\left|u\left(z,-\lambda^{2}\frac{1}{n}\right)\right|$
$\displaystyle\leq\frac{1}{\left(6\sqrt{n}+2-\varepsilon_{0}\right)^{2}}\left(\left|\left(6\sqrt{n}+1\right)^{2}\left(\lambda+\left|x_{0}\right|\right)\right|^{2}+R_{0}\right)$
$\displaystyle\leq\frac{1}{1+\varepsilon_{1}}\lambda^{2}\left[\left(1+\frac{\left|x_{0}\right|}{\lambda}\right)^{2}+\frac{R_{0}}{\lambda^{2}}\right],$

so, by choosing $\lambda>>\left|x_{0}\right|,R_{0}$, we obtain from (5.1)

$\sup_{x\in
B_{6\sqrt{n}+1}(0)}\left|\tilde{u}_{\lambda}(x,-\frac{1}{n})\right|\leq 1.$

We finally apply Proposition 4.2 to conclude

$D^{2}u(x_{0},0)=D^{2}u_{\lambda}(0,0)\leq C(n).$

Since this bound holds at any arbitrary point, it holds everywhere, and we may
apply [NY11, Proposition 2.1], which states that any convex ancient entire
solution with a uniform Hessian bound must be a quadratic polynomial. This
completes the proof.

## References

  * [BC03] Jiguang Bao and Jingyi Chen, _Optimal regularity for convex strong solutions of special Lagrangian equations in dimension 3_ , Indiana University mathematics journal (2003), 1231–1249.
  * [BCGJ03] Jiguang Bao, Jingyi Chen, Bo Guan, and Min Ji, _Liouville property and regularity of a Hessian quotient equation_ , American Journal of Mathematics 125 (2003), no. 2, 301–316.
  * [BDGM69] Enrico Bombieri, Ennio De Giorgi, and Mario Miranda, _An a priori bound related to minimal non-parametric hypersurfaces_ , Archive for Rational Mechanics and Analysis 32 (1969), no. 4, 255–267.
  * [Bha21] Arunima Bhattacharya, _Hessian estimates for Lagrangian mean curvature equation_ , Calculus of Variations and Partial Differential Equations 60 (2021), no. 6, 1–23.
  * [Bha22] by same author, _A note on the two-dimensional Lagrangian mean curvature equation_ , Pacific Journal of Mathematics 318 (2022), no. 1, 43–50.
  * [BMS24] Arunima Bhattacharya, Connor Mooney, and Ravi Shankar, _Gradient estimates for the Lagrangian mean curvature equation with critical and supercritical phase_ , arXiv:2205.13096, American Journal of Mathematics, to appear (2024).
  * [BW24] Arunima Bhattacharya and Jeremy Wall, _Hessian estimates for shrinkers, expanders, translators, and rotators of the Lagrangian mean curvature flow_ , arXiv:2403.07235, Preprint (2024).
  * [CCY12] Albert Chau, Jingyi Chen, and Yu Yuan, _Rigidity of entire self-shrinking solutions to curvature flows_ , Journal für die reine und angewandte Mathematik 2012 (2012), no. 664, 229–239.
  * [CCY13] by same author, _Lagrangian mean curvature flow for entire Lipschitz graphs II_ , Mathematische Annalen 357 (2013), no. 1, 165–183.
  * [CMI04] Tobias H. Colding and William P. Minicozzi II, _Sharp estimates for mean curvature flow of graphs_ , Journal für die reine und angewandte Mathematik 2004 (2004), no. 574, 187–195.
  * [CSY19] Jingyi Chen, Ravi Shankar, and Yu Yuan, _Regularity for convex viscosity solutions of special lagrangian equation_ , Communications on Pure and Applied Mathematics (2019).
  * [CWY09] Jingyi Chen, Micah Warren, and Yu Yuan, _A priori estimate for convex solutions to special Lagrangian equations and its application_ , Communications on Pure and Applied Mathematics 62 (2009), no. 4, 583–595.
  * [Din23] Qi Ding, _Liouville type theorems and Hessian estimates for special Lagrangian equations_ , Math. Ann. 386 (2023), no. 1-2, 1163–1200. MR 4585172
  * [EH91] Klaus Ecker and Gerhard Huisken, _Interior estimates for hypersurfaces moving by mean-curvarture_ , Inventiones mathematicae 105 (1991), no. 1, 547–569.
  * [ES92] Lawrence C Evans and Joel Spruck, _Motion of level sets by mean curvature III_ , The Journal of Geometric Analysis 2 (1992), no. 2, 121–150.
  * [HL82] Reese Harvey and H. Blaine Lawson, _Calibrated geometries_ , Acta Math. 148 (1982), 47–157.
  * [Kor86] Nicholas Korevaar, _An easy proof of the interior gradient bound for solutions to the prescribed mean-curvature equation_ , Proceedings of Symposia in Pure Mathematics, vol. 45, Amer Mathematical Soc 201 Charles st, Providence, RI 02940-2213, 1986, pp. 81–89.
  * [Li19] Caiyan Li, _A compactness approach to Hessian estimates for special Lagrangian equations with supercritical phase_ , Nonlinear Analysis 187 (2019), 434–437.
  * [LSS24] Jason D Lotay, Felix Schulze, and Gábor Székelyhidi, _Ancient solutions and translators of lagrangian mean curvature flow_ , Publications mathématiques de l’IHÉS (2024), 1–35.
  * [NY11] Tu A Nguyen and Yu Yuan, _A priori estimates for Lagrangian mean curvature flows_ , International Mathematics Research Notices 2011 (2011), no. 19, 4376–4383.
  * [Sha24] Ravi Shankar, _Hessian estimates for special lagrangian equation by doubling_ , arXiv preprint arXiv:2401.01034 (2024).
  * [SW02] Knut Smoczyk and Mu-Tao Wang, _Mean curvature flows of Lagrangian submanifolds with convex potentials_ , Journal of Differential Geometry 62 (2002), no. 2, 243–257.
  * [Tru72] Neil S. Trudinger, _A new proof of the interior gradient bound for the minimal surface equation in n dimensions_ , Proceedings of the National Academy of Sciences 69 (1972), no. 4, 821–823.
  * [TTW23] Chung-Jun Tsai, Mao-Pei Tsui, and Mu-Tao Wang, _Entire solutions of two-convex Lagrangian mean curvature flows_ , arXiv preprint arXiv:2309.09432 (2023).
  * [Wan02] Mu-Tao Wang, _Long-time existence and convergence of graphic mean curvature flow in arbitrary codimension_ , Inventiones mathematicae 148 (2002), no. 3, 525–543.
  * [Wan24] Jingbo Wan, _Sharp interior gradient estimate for area decreasing graphical mean curvature flow in arbitrary codimension_ , arXiv preprint arXiv:2407.05602 (2024).
  * [WY08] Micah Warren and Yu Yuan, _A liouville type theorem for special lagrangian equations with constraints_ , Communications in Partial Differential Equations 33 (2008), no. 5, 922–932.
  * [WY09a] by same author, _Explicit gradient estimates for minimal Lagrangian surfaces of dimension two_ , Math. Z. 262 (2009), no. 4, 867–879. MR 2511754
  * [WY09b] by same author, _Hessian estimates for the sigma-2 equation in dimension 3_ , Comm. Pure Appl. Math. 62 (2009), no. 3, 305–321. MR 2487850
  * [WY10] by same author, _Hessian and gradient estimates for three dimensional special Lagrangian equations with large phase_ , American Journal of Mathematics 132 (2010), no. 3, 751–770.
  * [WY14] Dake Wang and Yu Yuan, _Hessian estimates for special Lagrangian equations with critical and supercritical phases in general dimensions_ , American Journal of Mathematics 136 (2014), no. 2, 481–499.
  * [Zho22] Xingchen Zhou, _Hessian estimates to special Lagrangian equation on general phases with constraints_ , Calculus of Variations and Partial Differential Equations 61 (2022), no. 1, 1–14.


# Non-topological fractional fermion number in the Jackiw-Rossi model

Caio Almeida  Graduate Program in Physics (PPG-FIS), Universidade Federal do
ABC, Santo André, S.P., Brazil  Alberto Alonso-Izquierdo  Departamento de
Matematica Aplicada, University of Salamanca, Spain  Instituto de Fisica
Fundamental y Matematicas, IUFFyM, University of Salamanca, Spain  Rodrigo
Fresneda  CMCC-Universidade Federal do ABC, Santo André, S.P., Brazil  Juan
Mateos Guilarte  Departamento de Fisica Fundamental, University of Salamanca,
Spain  Instituto de Fisica Fundamental y Matematicas, IUFFyM, University of
Salamanca, Spain  Dmitri Vassilevich<EMAIL_ADDRESS>CMCC-Universidade
Federal do ABC, Santo André, S.P., Brazil  Physics Department, Tomsk State
University, Tomsk, Russia

###### Abstract

We compute the vacuum fermion current in $(2+1)$ dimensional Jackiw-Rossi
model by using the $1/m$ expansion. The current is expressed through a
weighted $\eta$-function with a matrix weight. In the presence of such a
weight, the usual proof of topological nature of $\eta(0)$ is not longer
applicable. Direct computations confirm the following surprising result: the
fermion number induced by vortices in the Jackiw-Rossi model is not
topological.

##  I Introduction

As we know from the pioneering paper by Jackiw and Rebbi [1] and from the
subsequent development reported in [2], the fermion number of solitons can
take fractional and even irrational values. In the known cases, the fermion
number is topological. This means that it depends on the boundary or
asymptotic values of the background fields and is not sensitive to smooth
variations of these fields in the interior of manifold. From the very
beginning, the fermion number fractionization had applications to condensed
matter physics [3]. More recently, this mechanism was applied to the physics
of topological insulators, see [4].

Among the planar ($2+1$-dimensional) solitonic systems, a prominent role is
played by the Abrikosov-Nielsen-Olesen (ANO) vortex. There are many possible
ways to couple fermions to this system, and thus there are many quantum
systems which include the ANO vortex as a bosonic sector. In a supersymmetric
model, the one-loop shift of the mass of the vortex was calculated in [5, 6].
In a pure bosonic model, this was done in [7, 8, 9] while (non-supersymmetric)
fermions were added in [10, 11].

The fermion number fractionization in $2+1$ dimensions on a pure gauge field
background was calculated in [12]. For a singular magnetic vortex this effect
was considered in [13]. For a pair of fermions coupled to both gauge and Higgs
fields of the ANO vertex, the half-integer fermion fractionization was
obtained in [14, 15] (see also the preceding papers [16, 17]). In these
models, the fermions have the elementary electric charge $e$ while the scalar
fields possess charge $2e$. Thus, the ANO vortex gets a fractional flux. There
is a way of coupling a single generation of fermions to the ANO system which
is given by the Jackiw-Rossi model [18]. This coupling reminds us of planar
superconducting systems. A candidate for the fractional flux vortex in such
systems was recently found experimentally [19]. This discovery motivated a
study [20] of fermion charge fractionization in the Jackiw-Rossi model. The
computations in this paper were based on the usual relation between vacuum
fermion number and the $\eta$-function of the Hamiltonian which is not correct
in the Jackiw-Rossi model, as will be demonstrated below.

The purpose of this paper is to analyze the vacuum fermion number in the
Jackiw-Rossi model paying special attention to its topological (or rather non-
topological) nature. First, we observe that the interaction between scalar and
spinor fields does not allow to immediately relate the fermion number to the
$\eta$-function of an operator. One has to double the spinor components. This
is similar to what has been done in [18, 21] to analyses the zero-energy
spectrum, but we do this in the path integral formalism following the method
of [22]. The fermion density is then related to an $\eta$-function which,
however, is weighted with a matrix. The presence of this matrix destroys the
standard proof of vanishing local variations of $\eta$ and of the topological
nature of this quantity. We go on by computing a few leading term in the large
mass expansion of the fermion number with the heat kernel methods and confirm
the presence of non-topological contributions depending on the profiles of
magnetic and Higgs fields rather than on their global characteristics.

This paper is organized as follows. In the next section, we derive an
expression for the vacuum fermion number in Jackiw-Rossi model in terms of an
$\eta$ function with a matrix weight. Section III is dedicated to the heat
kernel evaluation of fermion number. First, we show why the standard proof of
topological nature does not work for weighted $\eta$ function. Then, we
pinpoint the non-topological fermion number within the large mass expansion.
Some concluding remarks are given in the last section.

##  II Fermion number in the Jackiw-Rossi model

The Jackiw-Rossi model [18] in $(2+1)$ dimensions is described by the
Lagrangian

$\mathcal{L}=\bar{\psi}(\gamma^{\mu}(\mathrm{i}\partial_{\mu}-eA_{\mu}))\psi-\frac{1}{2}\mathrm{i}g\phi\bar{\psi}\psi^{C}+\frac{1}{2}\mathrm{i}g^{*}\phi^{*}\bar{\psi}^{C}\psi-m\bar{\psi}\psi$
(1)

governing the dynamics of a two complex component spinor field $\psi$ coupled
to a gauge and a complex scalar fields $A_{\mu}$, $\phi$.

As compared to the original work [18] a mass term has been added. For
convenience, let us take the $\gamma$-matrices in Majorana representation:

$\gamma^{0}=\sigma^{2},\qquad\gamma^{1}=\mathrm{i}\sigma^{1},\qquad\gamma^{2}=\mathrm{i}\sigma^{3}.$
(2)

Then the charge conjugation matrix can be taken as $C=-\gamma^{0}$. We have
the usual relations $C\gamma^{\mu}C^{-1}=-\gamma^{\mu\,T}$,
$\psi^{C}=\psi^{*}$, etc.

We assume that bosonic fields $A_{\mu}$ and $\phi$ belong to the topological
class of an ANO vortex. This configuration is static, so that $A_{0}=0$ and
all fields do not depend on time. We are not going to use the exact profile
functions, though it will be important to us that this configuration is
localized somewhere near the origin. If $r$ is the radial coordinate, for
$r\to\infty$ we have

$|\phi|\to v,\qquad D_{j}\phi\to 0,\qquad F_{jk}\to 0.$ (3)

Here and in what follows $x^{j}$, $x^{k}$, etc denote spatial coordinates.
$D_{j}\phi=(\partial_{j}+2\mathrm{i}eA_{j})\phi$ is a gauge covariant
derivative, depending on the charge of the field it acts upon; therefore, in
our notation, $D_{j}\phi^{*}\equiv(D_{j}\phi)^{*}$. Note that the electric
charge of $\phi$ is $2e$, $F_{jk}\equiv\partial_{j}A_{k}-\partial_{k}A_{j}$,
and $v$ is a minimum of the Higgs potential. All functions in (3) go to their
asymptotic values exponentially fast. Let $\mathcal{N}\in\mathbb{Z}$ be the
topological charge of the vortex. The magnetic flux quantization condition

$\frac{e}{\pi}\int d^{2}x\,F_{12}=\mathcal{N}$ (4)

has an unusual factor on the right hand side due to the charge $2e$ of $\phi$.
This is why we say that the vortex has a fractional flux.

The Lagrangian (1) besides the $\bar{\psi}\psi$ contains also the
$\bar{\psi}\psi^{*}$ and $\bar{\psi}^{*}\psi$ couplings to the Higgs field.
Thus, it does not have the form that allows to relate immediately the states
to eigenfunctions of some differential operator. To overcome this difficulty
we pass to doubled spinors following the approach developed in the paper [22]
(see also [23, 24]). We introduce

$\Psi:=\begin{pmatrix}\psi\\\ \bar{\psi}^{T}\end{pmatrix}.$ (5)

With the help of identities

$\displaystyle\int
d^{3}x\bar{\psi}\gamma^{\mu}(\mathrm{i}\partial_{\mu}-eA_{\mu})\psi=\int
d^{3}x\psi^{T}\gamma^{\mu\,T}(\mathrm{i}\partial_{\mu}+eA_{\mu})\bar{\psi}^{T},$
$\displaystyle\bar{\psi}^{C}\psi=\psi^{T}\gamma^{0}\psi,\qquad\bar{\psi}\psi^{C}=-\bar{\psi}\gamma^{0}\bar{\psi}^{T},\qquad-m\bar{\psi}\psi=m\psi^{T}\bar{\psi}^{T}$

we rewrite the action as

$S=\frac{1}{2}\int d^{3}x\Psi^{T}\widehat{F}\Psi$ (6)

with

$\widehat{F}=\begin{pmatrix}\mathrm{i}g^{*}\phi^{*}\gamma^{0}&\gamma^{\mu\,T}(\mathrm{i}\partial_{\mu}+eA_{\mu})+m\\\
\gamma^{\mu}(\mathrm{i}\partial_{\mu}-eA_{\mu})-m&\mathrm{i}g\phi\gamma^{0}\end{pmatrix}.$
(7)

The corresponding Hamiltonian reads

$H=\begin{pmatrix}\alpha^{j}(\mathrm{i}\partial_{j}-eA_{j})-\beta
m&\mathrm{i}g\phi\\\
-\mathrm{i}g^{*}\phi^{*}&-\alpha^{j}(\mathrm{i}\partial_{j}+eA_{j})-\beta
m\end{pmatrix}.$ (8)

Here, as usual, $\beta\equiv\gamma^{0}$ and $\alpha^{j}=\beta\gamma^{j}$.

Let us consider the effective action $W$ which is obtained by integrating out
the fermionic degrees of freedom,

$e^{\mathrm{i}W}=\int\mathcal{D}\psi\,\mathcal{D}\bar{\psi}\,\exp\left(\mathrm{i}\int
d^{3}x\mathcal{L}\right).$ (9)

This action depends on the background bosonic fields $\phi$ and $A_{\mu}$. The
charge density is given by the variational derivative

$j^{0}=-\frac{1}{e}\,\frac{\delta W}{\delta A_{0}}.$ (10)

The same effective action $W$ can be written through a path integral over the
doubled spinors $\Psi$ as

$W=-\mathrm{i}\,\ln\int\mathcal{D}\Psi\,\exp\left(\mathrm{i}\tfrac{1}{2}\int
d^{3}x\,\Psi^{T}\widehat{F}\Psi\right)=-\frac{\mathrm{i}}{2}\ln\det(\widehat{F}),$
(11)

see [22]. The functional integration measure became
$\mathcal{D}\Psi=\mathcal{D}\psi\,\mathcal{D}\bar{\psi}$.

Symbolically, we may write

$j^{0}=\frac{\mathrm{i}}{2e}\,\mathrm{Tr}\left(\frac{\delta\widehat{F}}{\delta
A_{0}}\,\widehat{F}^{-1}\right).$ (12)

To give precise meaning to this formula one has to invert $\widehat{F}$ and
regularize the functional trace.

After having calculated the variational derivative in (12) one puts the
background fields to their values for the static vortex configuration. On such
a background, the eigenfunctions of Hamiltonian (8) can be taken depending on
the spatial coordinates $\vec{x}$ only,

$H\Psi_{n}(\vec{x})=E_{n}\Psi_{n}(\vec{x}).$ (13)

The energy spectrum has both discrete and continuous parts. To avoid notation
clutter we write the formulas below as if the whole spectrum were discrete.

For $g=0$, the Hamiltonian $H$ consists of two hermitian anticommuting parts.
Thus one can easily show that $(H(g=0))^{2}\geq m^{2}$. Consequently, if
$|m|>|g\phi|$ the full Hamiltonian does not have zero energy eigenstates.

The vectors

$\Psi_{\omega,n}(\vec{x},t)=(2\pi)^{-1/2}e^{-\mathrm{i}\omega
t}\,\Psi_{n}(\vec{x})$ (14)

form a basis for the space of square integrable 4-spinors on $\mathbb{R}^{3}$.
To compute (12) one has to sandwich the expression under the trace between
$\Psi_{\omega,n}^{\dagger}$ and $\Psi_{\omega,n}$, integrate over $\omega$ and
sum over $n$. To regularize the $\omega$-integral we use a symmetric time-
splitting regularization. Namely, we take two eigenvectors at shifted time
arguments, $\Psi_{\omega,n}^{\dagger}(\vec{x},t)$ and
$\Psi_{\omega,n}(\vec{x},t+\Delta t)$. After computing the integral, we take
the limits $\tfrac{1}{2}(\lim_{\Delta t\to+0}+\lim_{\Delta t\to-0})$. To
regularize the sum, we multiply the expression by $|E_{n}|^{-s}$ with $\Re s$
sufficiently large to ensure the convergence and analytically continue to
$s=0$ afterwards.

$\displaystyle
j^{0}(\vec{x},t)=-\frac{\mathrm{i}}{2}\frac{1}{2}\left(\lim_{\Delta
t\to+0}+\lim_{\Delta
t\to-0}\right)\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\sum_{n}|E_{n}|^{-s}$
$\displaystyle\qquad\qquad\times\Psi_{n}^{\dagger}(\vec{x})\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\Psi_{n}(\vec{x})\frac{e^{-\mathrm{i}\omega\,\Delta
t}}{\omega+\mathrm{i}\,0\,\mathrm{sgn}(\omega)-E_{n}}\,.$ (15)

After performing the integration over $\omega$ one obtains

$j^{0}(\vec{x},t)=-\frac{1}{4}\sum_{n}\Psi_{n}^{\dagger}(\vec{x})\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\Psi_{n}(\vec{x})\,\mathrm{sgn}\,(E_{n})\,|E_{n}|^{-s}$ (16)

The analytic continuation to $s=0$ is understood in both formulas (15) and
(16).

Let $Q$ be a smooth bounded matrix-valued function (a smooth endomorphism).
The $\eta$ function of $H$ smeared with $Q$ is defined as

$\eta(s,H;Q)=\mathrm{Tr}\,\left(Q\,\mathrm{sgn}\,(H)\,|H|^{-s}\right)=\mathrm{Tr}\,\left(Q\cdot(H^{2})^{-s/2}\cdot(H/|H|)\right)=\mathrm{Tr}\,\left(Q\cdot(H^{2})^{-\frac{s+1}{2}}H\right).$
(17)

Here again $s$ is a complex parameter. The trace in (17) exists if $\Re\,s$ is
sufficiently large. This function can be analytically continued as a
meromorphic function to the whole complex plane. At $s=0$, equation (16)
yields

$\int d^{2}x\,j^{0}(x)\,\rho(x)=-\frac{1}{4}\,\eta(0,H;\rho\tau_{3})\,,$ (18)

where $\rho$ is a smooth localizing function of compact support, and
$\tau_{3}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right)$. An
integrated version of (16) gives an expression for the fermion number $N$ in
terms of the $\eta$ function,

$N\equiv\int d^{2}x\,j^{0}(x)=-\frac{1}{4}\,\eta(0,H;\tau_{3}).$ (19)

There are two important differences from the corresponding formula derived in
the seminal paper [12]. These are the coefficient $1/4$ instead of $1/2$ and
the presence of $\tau_{3}$ in the $\eta$ function. Both are caused by our
spinor field doubling procedure. The presence of $\tau_{3}$ has a profound
consequence: the standard proof that $\eta(0)$ is topological in 2D does not
work any more.

##  III Heat kernel computations of the fermion current

###  III.1 Why the standard proof of $N$ being topological does not work for
the JR model

Here we study local variations of the $\eta$-function with and without a
matrix weighting factor. Our method goes back to the paper by Atiyah, Patodi
and Singer[25]. We closely follow the procedure presented in [26]. A slightly
different method was used in [27].

Let $H(\varepsilon)=H+\varepsilon\,h$ where $h$ is a perturbation caused by an
infinitesimal localized variation of background bosonic fields $\phi$ and $A$.
Let us consider the case $Q=1$. By using Lemma 1.10.2 of [26] we can express
the variation of the $\eta$ function through a $\zeta$ function weighted with
$h$

$\left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\right|_{\varepsilon=0}\eta(s,H(\varepsilon))=-s\mathrm{Tr}\,\left[h\,(H^{2})^{-\frac{s+1}{2}}\right]=-s\zeta\left(\frac{s+1}{2},H^{2};h\right).$
(20)

Now, we need some basics on the spectral functions. Let $L$ be a Laplace type
operator on a manifold $M$ of dimension $n$ with or without boundary. Let $h$
be a smooth matrix valued function. Then, residues of the $\zeta$ function can
be expressed by the formula

$\mathrm{Res}_{u=\frac{n-k}{2}}\left(\Gamma(u)\zeta(u,L;h)\right)=a_{k}(L;h)$
(21)

through the heat kernel coefficients defined through the following asymptotic
expansion at $t\to+0$

$\mathrm{Tr}\,\left(he^{-tL}\right)\simeq\sum_{k=0}^{\infty}t^{\frac{k-n}{2}}a_{k}(L;h).$
(22)

By Eq. (20), the derivative
$(\mathrm{d}\eta(0,H))/(\mathrm{d}\varepsilon)|_{\varepsilon=0}$ is given by
the residue of $\zeta(u,H^{2};h)$ function at $u=\tfrac{1}{2}$ which is in
turn proportional to $a_{1}(H^{2};h)$. Since $h$ is localized inside the
manifold and does not extend to boundaries or asymptotic regions, the
coefficient $a_{1}(H^{2};h)$ vanishes. We conclude that $\eta(0,H)$ does not
change under local variations $H\to H(\varepsilon)=H+\varepsilon h$ and thus
is a topological invariant.

The key point of the proof presented above was the Eq. (20) relating the
variation of the $\eta$ function to a residue of the $\zeta$ function which
happened to be local and vanishing in the dimension $n=2$. Roughly speaking,
to get (20) one needs to differentiate $\eta(0,H)$ as if it were a usual
function of a commutative variable. This property is ensured by the
possibility of reordering operators under the trace. This possibility is
(partially) lost if $Q$ does not commute with $H$. In such a case, the
variation of $\eta(0,H;Q)$ cannot be written in the simple form of (20) and
all subsequent arguments break down.

###  III.2 Computations for the JR model

To evaluate the large mass expansion of the current (18) we shall use the
method proposed in [28, 29]. With the help of the identity

$\int_{0}^{\infty}dt\,t^{a}e^{-bt}=b^{-(1+a)}\Gamma(1+a)$ (23)

we write

$\eta(s,H,\rho\tau_{3})=\frac{1}{\Gamma\left(\frac{s+1}{2}\right)}\int_{0}^{\infty}dt\,t^{\frac{s-1}{2}}\mathrm{Tr}\,\left(\rho\tau_{3}He^{-tH^{2}}\right).$
(24)

Let us introduce a shifted operator $H_{\rho}=H-\varepsilon\rho\tau_{3}$ with
$\varepsilon$ being a real parameter. Then

$\eta(s,H,\rho\tau_{3})=\frac{1}{2\Gamma\left(\frac{s+1}{2}\right)}\int_{0}^{\infty}dt\,t^{\frac{s-3}{2}}\left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\right|_{\varepsilon=0}\mathrm{Tr}\,\left(e^{-tH_{\rho}^{2}}\right).$
(25)

To evaluate this expression by using a large mass expansion we isolate $m^{2}$
in $H_{\rho}^{2}$ and take the limit $s\to 0$ to obtain

$\eta(0,H,\rho\tau_{3})=\frac{1}{2\sqrt{\pi}}\int_{0}^{\infty}dt\,t^{-\frac{3}{2}}\left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\right|_{\varepsilon=0}e^{-tm^{2}}\mathrm{Tr}\,\left(e^{-t\tilde{H}_{\rho}^{2}}\right),$
(26)

where $\tilde{H}_{\rho}^{2}\equiv H_{\rho}^{2}-m^{2}$. Next, we make the heat
kernel expansion (22) and integrate over $t$.

$\displaystyle\eta(0,H;\rho\tau_{3})$ $\displaystyle\simeq$
$\displaystyle\frac{1}{2\sqrt{\pi}}\int_{0}^{\infty}dt\,\sum_{k=0}^{\infty}t^{\frac{k-5}{2}}\left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\right|_{\varepsilon=0}a_{k}(\tilde{H}_{\rho}^{2})\,e^{-tm^{2}}$
(27) $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{\pi}}\sum_{k}\Gamma\left(\frac{k-3}{2}\right)\,|m|^{3-k}\left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\right|_{\varepsilon=0}a_{k}(\tilde{H}_{\rho}^{2}).$

Here $a_{k}(\tilde{H}_{\rho}^{2})\equiv a_{k}(\tilde{H}_{\rho}^{2};1)$. The
integral above is convergent if the contributions of heat kernel coefficients
$a_{k}$ with $k\leq 3$ vanish. We shall check this condition below.

To be able to use universal expressions for the heat kernel coefficients (see,
e.g., [30]) we represent the operator $\tilde{H}_{\rho}^{2}$ in the canonical
form

$\tilde{H}_{\rho}^{2}=-(\nabla_{j}\nabla_{j}+E),$ (28)

where $\nabla_{j}=\partial_{j}+\omega_{j}$ plays the role of a covariant
derivative while $E$ is a matrix valued potential. For our operator they read

$\displaystyle
E=\begin{pmatrix}\frac{e}{2}\beta\,\epsilon^{jk}F_{jk}-|g\phi|^{2}&g\alpha^{j}D_{j}\phi+2\mathrm{i}\beta
g\phi m\\\ g^{*}\alpha^{j}D_{j}\phi^{*}-2\mathrm{i}\beta
g^{*}\phi^{*}m&-\frac{e}{2}\beta\,\epsilon^{jk}F_{jk}-|g\phi|^{2}\end{pmatrix}-2\varepsilon\rho\beta
m\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$ (29)
$\displaystyle\omega_{j}=\begin{pmatrix}\mathrm{i}eA_{j}&0\\\
0&-\mathrm{i}eA_{j}\end{pmatrix}+\mathrm{i}\alpha^{j}\varepsilon\rho$ (30)

In this section, we are working in a Euclidean space with a positive unit
metric. We still keep the distinction between upper and lower indices of some
quantities which have a $(2+1)$-dimensional origin. For example, $A$ always
appears with a subscript, while $\alpha$ and $\gamma$ come with superscripts.
The summation over repeated indices is always done with the Kronecker symbol
independently of the position of indices. This prescription destroys the
balance between upper and lower indices within formulas, but keeps the
notations simple and unambiguous.

Each heat kernel coefficient $a_{k}$ is an integral of a trace of a local
polynomial constructed from $E$, the field strength
$\Omega_{ij}=[\nabla_{i},\nabla_{j}]$, and their repeated covariant
derivatives (e.g., $E_{;j}=[\nabla_{j},E]$, etc). For example,

$\displaystyle
E_{;j}=\begin{pmatrix}\frac{e}{2}\beta\partial_{j}(\epsilon^{kl}F_{kl})-2\varepsilon\beta
m\partial_{j}\rho-\partial_{j}|g\phi|^{2}&g\alpha^{k}(D_{j}D_{k}\phi)+2\mathrm{i}g\beta
m(D_{j}\phi)\\\ g^{*}\alpha^{k}(D_{j}D_{k}\phi^{*})-2\mathrm{i}g^{*}\beta
m(D_{j}\phi^{*})&-\frac{e}{2}\beta\partial_{j}(\epsilon^{kl}F_{kl})+2\varepsilon\beta
m\partial_{j}\rho-\partial_{j}|g\phi|^{2}\end{pmatrix}$
$\displaystyle\qquad+\begin{pmatrix}-\mathrm{i}\gamma^{j}e\hskip
1.5pt\epsilon^{kl}F_{kl}&2g\beta\epsilon^{jk}(D_{k}\phi)+4\gamma^{j}mg\phi\\\
2g^{*}\beta\epsilon^{jk}(D_{k}\phi^{*})-4\gamma^{j}mg^{*}\phi^{*}&\mathrm{i}\gamma^{j}e\hskip
1.5pt\epsilon^{kl}F_{kl}\end{pmatrix}\varepsilon\rho$ (31)

and

$\Omega_{ij}=\begin{pmatrix}\mathrm{i}eF_{ij}+\mathrm{i}\varepsilon\left(\alpha^{j}(\partial_{i}\rho)-\alpha^{i}(\partial_{j}\rho)\right)&0\\\
0&-\mathrm{i}eF_{ij}+\mathrm{i}\varepsilon\left(\alpha^{j}(\partial_{i}\rho)-\alpha^{i}(\partial_{j}\rho)\right)\end{pmatrix}\,.$
(32)

All invariants entering $a_{k}$ have the canonical mass dimension $k$. On
manifolds without boundaries, all coefficients with odd values of $k$ vanish.

By a direct computation with the expressions from [30], one obtains that the
contributions of $a_{0}$ and $a_{2}$ to (27) vanish thus fulfilling the
consistency condition presented below (27).

Let us compute the current as an expansion in $\phi$ and its derivatives
keeping the terms up to $D^{2}$ and $\phi^{2}$. Since $[D_{j},D_{k}]\propto
F_{jk}$ we shall also keep the terms with $F$ and $F\phi^{2}$, while $m$ can
enter with any power.

It is important to establish upper bounds on the number $k$ of the heat kernel
coefficient which contains the required invariants. Consider the term
$\rho\epsilon^{jk}F_{jk}$. It has canonical mass dimension $3$. Thus, in
$a_{k}$ it has to be multiplied by $m^{k-3}$. Since $\omega$ does not contain
$m$, this requires product of $E$ or of its derivatives at least $k-3$ times –
an expression which has the mass dimension greater than or equal to $2(k-3)$.
Since the mass dimension of $a_{k}$ is $k$, we have the upper bound $k\leq 6$.
In a similar way one comes to the conclusion that $\rho(D\phi)(D\phi^{*})$ and
$\rho|g\phi|^{2}F$ terms may appear for $k\leq 10$. By refining these
arguments one can exclude a lot of possible terms in the expansion and even
improve the bounds mentioned above. At any rate, with explicit expressions
from [31] for flat space heat kernel coefficients up to $a_{12}$, the rest may
be done by a Wolfram Mathematica script.

We obtain that, up to the order considered, just a few terms in the heat
kernel expansion contribute. The result reads

$\displaystyle a_{4}=\frac{1}{4\pi}\int
d^{2}x\,\mathrm{tr}\left(\tfrac{1}{2}E^{2}+\dots\right)=-\frac{1}{\pi}\int
d^{2}x\,e\epsilon^{jk}F_{jk}\varepsilon\rho m+\dots,$ (33) $\displaystyle
a_{6}=\frac{1}{4\pi}\int
d^{2}x\,\mathrm{tr}\left(\tfrac{1}{6}E^{3}-\tfrac{1}{12}E_{;j}E_{;j}+\dots\right)=\frac{1}{\pi}\int
d^{2}x\,\varepsilon\rho
m\,|g|^{2}\epsilon^{jk}\left(\frac{\mathrm{i}}{3}(D_{j}\phi)(D_{k}\phi^{*})+\frac{5e}{3}|\phi|^{2}F_{jk}\right)+\dots,$
(34) $\displaystyle a_{8}=\frac{1}{4\pi}\int
d^{2}x\,\mathrm{tr}\left(\tfrac{1}{24}E^{4}+\dots\right)=-\frac{2}{3\pi}\int
d^{2}x\,em^{3}|g\phi|^{2}\varepsilon\rho\,\epsilon^{jk}F_{jk}+\dots\,,$ (35)

where dots denote irrelevant terms.

Thus, in the approximation adopted here,

$j^{0}=\frac{1}{8\pi}\,\frac{m}{|m|}\left[e\epsilon^{jk}F_{jk}-\frac{\mathrm{i}|g|^{2}}{6m^{2}}\,\epsilon^{jk}(D_{j}\phi)(D_{k}\phi^{*})-\frac{|g\phi|^{2}}{3m^{2}}e\epsilon^{jk}F_{jk}\right].$
(36)

The integral of $j^{0}$ gives the vacuum fermion number

$N=\frac{\mathcal{N}}{4}\,\frac{m}{|m|}-\frac{|g|^{2}e}{48\pi\,m\,|m|}\int
d^{2}x\,|\phi|^{2}\epsilon^{jk}F_{jk}.$ (37)

To obtain this expression we integrated by parts and used the asymptotic
conditions (3) together with the relation (4). The first term on right hand
side of (37) describes the (expected) quarter-integer quantization of the
fermion number in the absence of scalar field $\phi$. The second term depends
on the profiles of $|\phi|^{2}$ and $F_{jk}$ in the interior of manifold and
thus is not topological.

##  IV Conclusions

In this paper, we have expressed the vacuum fermion number of the Jackiw-Rossi
model through an $\eta$-invariant of a matrix-weighted Hamiltonian. We have
pinpointed the reason why the standard proof of the topological nature of the
(fractional) fermion number fails and we have also explicitly computed a non-
topological contribution to this quantity.

We have computed the fermion current in just a few leading orders of the large
mass expansion, since this was enough for our purposes. If needed, further
terms can also be calculated with the help of the flat space heat kernel
expansion from the paper [31]. One can also use resummations of the heat
kernel, see [32, 33].

The fact that the fermion number depends on the profiles of the magnetic field
and of the Higgs field should have some consequences for condensed matter
physics. We are not ready to go deeper into this subject. We just mention a
potentially related work which studies, both theoretically and experimentally,
the influence of non-uniformity of the magnetic field on Hall conductivity for
various planar systems [34].

Speaking about future prospects, we would also like to mention the work [35]
which studies relations between the parameters of solitons and the fermion
spectrum. Probably, these results can be lifted to the quantum level to gain
information about the fermion fractionization and other similar effects.

###### Acknowledgements.

This work was done in the framework of an agreement between the São Paulo
Research Foundation (FAPESP) and the University of Salamanca, project
2017/50294-1 (SPRINT). The work was supported in parts by the project
2016/03319-6 of FAPESP, by the grants 305594/2019-2 and 428951/2018-0 of CNPq.
Besides, D.V. was supported by the Tomsk State University Competitiveness
Improvement Program. AAI and JMG acknowledge the Junta de Castilla y Leon for
partial financial support, under Grants No. BU229P18 and No. SA067G19. CA
received a graduate scholarship directly provided by the Federal University of
ABC and, hence, would like to thank the institution for the financial support.

## References

  * [1] R. Jackiw and C. Rebbi, “Solitons with Fermion Number 1/2,” Phys. Rev. D 13, 3398 (1976). doi:10.1103/PhysRevD.13.3398
  * [2] A. J. Niemi and G. W. Semenoff, “Fermion Number Fractionization in Quantum Field Theory,” Phys. Rept. 135, 99 (1986). doi:10.1016/0370-1573(86)90167-5
  * [3] R. Jackiw and J. R. Schrieffer, “Solitons with Fermion Number 1/2 in Condensed Matter and Relativistic Field Theories,” Nucl. Phys. B 190, 253 (1981). doi:10.1016/0550-3213(81)90557-5
  * [4] X. L. Qi and S. C. Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, no. 4, 1057 (2011) doi:10.1103/RevModPhys.83.1057 [arXiv:1008.2026 [cond-mat.mes-hall]].
  * [5] D. V. Vassilevich, “Quantum corrections to the mass of the supersymmetric vortex,” Phys. Rev. D 68, 045005 (2003) doi:10.1103/PhysRevD.68.045005 [hep-th/0304267].
  * [6] A. Rebhan, P. van Nieuwenhuizen and R. Wimmer, “Nonvanishing quantum corrections to the mass and central charge of the N = 2 vortex and BPS saturation,” Nucl. Phys. B 679, 382 (2004) doi:10.1016/j.nuclphysb.2003.11.030 [hep-th/0307282].
  * [7] A. Alonso Izquierdo, W. Garcia Fuertes, M. de la Torre Mayado and J. Mateos Guilarte, “Quantum corrections to the mass of self-dual vortices,” Phys. Rev. D 70, 061702 (2004) doi:10.1103/PhysRevD.70.061702 [hep-th/0406129].
  * [8] A.Alonso Izquierdo, W. Garcia Fuertes, J. Mateos Guilarte, M. de la Torre Mayado, “Quantum oscillations of self-dual Abrikosov-Nielsen-Olesen vortices”, Phys. Rev. D 71, 125010 [arXiv:hep-th/0504143].
  * [9] A.Alonso Izquierdo, J. Mateos Guilarte, M. de la Torre Mayado, “Quantum magnetic flux lines, BPS vortex zero modes, and one-loop string tension shifts”, Phys. Rev. D 94, 045008 (2016) [arXiv:1605.09175[hep-th]].
  * [10] M.Bordag and I. Drozdov, “Fermionic vacuum energyfrom a Nielsen-Olesen vortex”, Phys Rev. D 68, 065026 (2003)[arXiv:hep-th/0305002]
  * [11] N. Graham, V. Khemani, M. Quandt, O. Schroeder and H. Weigel, “Quantum QED flux tubes in 2+1 and 3+1 dimensions,” Nucl. Phys. B 707, 233 (2005) doi:10.1016/j.nuclphysb.2004.11.057 [hep-th/0410171].
  * [12] A. J. Niemi and G. W. Semenoff, “Axial Anomaly Induced Fermion Fractionization and Effective Gauge Theory Actions in Odd Dimensional Space-Times,” Phys. Rev. Lett. 51, 2077 (1983). doi:10.1103/PhysRevLett.51.2077
  * [13] Y. A. Sitenko, “Selfadjointness of the Dirac Hamiltonian and fermion number fractionization in the background of a singular magnetic vortex,” Phys. Lett. B 387, 334 (1996) doi:10.1016/0370-2693(96)01042-8 [hep-th/9604143].
  * [14] C. Chamon, C. Y. Hou, R. Jackiw, C. Mudry, S. Y. Pi and A. P. Schnyder, “Irrational vs. rational charge and statistics in two-dimensional quantum systems,” Phys. Rev. Lett. 100, 110405 (2008) doi:10.1103/PhysRevLett.100.110405 [arXiv:0707.0293 [cond-mat.str-el]].
  * [15] C. Chamon, C. Y. Hou, R. Jackiw, C. Mudry, S. Y. Pi and G. Semenoff, “Electron fractionalization for two-dimensional Dirac fermions,” Phys. Rev. B 77, 235431 (2008) doi:10.1103/PhysRevB.77.235431 [arXiv:0712.2439 [hep-th]].
  * [16] C. Y. Hou, C. Chamon and C. Mudry, “Electron fractionalization in two-dimensional graphene-like structures,” Phys. Rev. Lett. 98, 186809 (2007) doi:10.1103/PhysRevLett.98.186809 [cond-mat/0609740 [cond-mat.mes-hall]].
  * [17] R. Jackiw and S.-Y. Pi, “Chiral gauge theory for graphene,” Phys. Rev. Lett. 98, 266402 (2007) doi:10.1103/PhysRevLett.98.266402 [cond-mat/0701760 [cond-mat.str-el]].
  * [18] R. Jackiw and P. Rossi, “Zero Modes of the Vortex - Fermion System,” Nucl. Phys. B 190, 681 (1981). doi:10.1016/0550-3213(81)90044-4
  * [19] Y. Tanaka, H. Yamamori, T. Yanagisawa, T. Nishio and S. Arisawa, ”Experimental formation of a fractional vortex in a superconducting bi-layer,” Physica C: Superconductivity and its Applications 548, 44 (2018). doi:10.1016/j.physc.2018.02.001
  * [20] T. Yanagisawa, ”Zero-Energy Modes, Fractional Fermion Numbers and The Index Theorem in a Vortex-Dirac Fermion System,” Symmetry 12, 373 (2020). doi:10.3390/sym12030373
  * [21] E. J. Weinberg, “Index Calculations for the Fermion-Vortex System,” Phys. Rev. D 24, 2669 (1981). doi:10.1103/PhysRevD.24.2669
  * [22] Y. Novozhilov, A. Pronko and D. Vassilevich, “Quark path integral with diquark currents and its anomalies,” Phys. Lett. B 321, 425 (1994). doi:10.1016/0370-2693(94)90270-4
  * [23] R. D. Ball, “Majorana Fermion Determinants,” Phys. Lett. B 227, 445 (1989). doi:10.1016/0370-2693(89)90959-3
  * [24] Y. Novozhilov, A. Pronko and D. Vassilevich, “Extended chiral transformations including diquark fields as parameters,” Phys. Lett. B 343, 358 (1995) Erratum: [Phys. Lett. B 351, 601 (1995)] doi:10.1016/0370-2693(95)00391-W, 10.1016/0370-2693(94)01451-H [hep-ph/9411421].
  * [25] M. F. Atiyah, V. K. Patodi and I. M. Singer, “Spectral asymmetry and Riemannian geometry. III,” Math. Proc. Cambridge Phil. Soc. 79, 71 (1976). doi:10.1017/S0305004100052105
  * [26] P. B. Gilkey, Invariance theory, the heat equation, and the Atiyah-Singer index theorem (Publish or Perish, Wilmington, 1984).
  * [27] L. Alvarez-Gaume, S. Della Pietra and G. W. Moore, “Anomalies and Odd Dimensions,” Annals Phys. 163, 288 (1985). doi:10.1016/0003-4916(85)90383-5
  * [28] A. Alonso-Izquierdo, R. Fresneda, J. Mateos Guilarte and D. Vassilevich, “Soliton Fermionic number from the heat kernel expansion,” Eur. Phys. J. C 79, no. 6, 525 (2019) doi:10.1140/epjc/s10052-019-7041-8 [arXiv:1905.09030 [hep-th]].
  * [29] J. Mateos Guilarte and D. Vassilevich, “Fractional fermion number and Hall conductivity of domain walls,” Phys. Lett. B 797, 134935 (2019) doi:10.1016/j.physletb.2019.134935 [arXiv:1908.07989 [hep-th]].
  * [30] D. V. Vassilevich, “Heat kernel expansion: User’s manual,” Phys. Rept. 388, 279 (2003) doi:10.1016/j.physrep.2003.09.002 [hep-th/0306138].
  * [31] D. Fliegner, P. Haberl, M. G. Schmidt and C. Schubert, “The Higher derivative expansion of the effective action by the string inspired method. Part 2,” Annals Phys. 264, 51 (1998) doi:10.1006/aphy.1997.5778 [hep-th/9707189].
  * [32] A. O. Barvinsky and G. A. Vilkovisky, “Covariant perturbation theory. 2: Second order in the curvature. General algorithms,” Nucl. Phys. B 333, 471 (1990). doi:10.1016/0550-3213(90)90047-H
  * [33] I. G. Avramidi, “Covariant techniques for computation of the heat kernel,” Rev. Math. Phys. 11, 947 (1999) doi:10.1142/S0129055X99000295 [hep-th/9704166].
  * [34] J. Schirmer, R. Kumar, V. Bagwe, P. Raychaudhuri, T. Taniguchi, K. Watanabe, C.-X. Liu, A. Das and J. K. Jain, ”Hall effect for Dirac electrons in graphene exposed to an Abrikosov flux lattice,” EPL (Europhysics Letters) 132, 37002 (2020). doi:10.1209/0295-5075/132/37002 [arXiv:2010.14400]
  * [35] D. Bazeia, A. Mohammadi and D. C. Moreira, “Fermions in the presence of topological structures under geometric constrictions,” Phys. Rev. D 103, no. 2, 025003 (2021) doi:10.1103/PhysRevD.103.025003 [arXiv:2009.00737 [hep-th]].