# Gluing invariants of Donaldson–Thomas type
Part I: the Darboux stack

###### Abstract.

Let $X$ be a $(-1)$–shifted symplectic derived Deligne–Mumford stack. In this
paper we introduce the Darboux stack of $X$, parametrizing local presentations
of $X$ as a derived critical locus of a function $f$ on a smooth formal scheme
$U$. Local invariants such as the Milnor number $\mu_{f}$, the perverse sheaf
of vanishing cycles $\PJoyce_{U,f}$ and the category of matrix factorizations
$\MF(U,f)$ are naturally defined on the Darboux stack, without ambiguity. The
stack of non-degenerate flat quadratic bundles acts on the Darboux stack and
our main theorem is the contractibility of the quotient stack when taking a
further homotopy quotient identifying isotopic automorphisms. As a corollary
we recover the gluing results for vanishing cycles by
Brav–Bussi–Dupont–Joyce–Szendrői.

In a second part (to appear), we will apply this general mechanism to glue the
motives of the locally defined categories of matrix factorizations $\MF(U,f)$
under the prescription of additional orientation data, thus answering
positively conjectures by Kontsevich–Soibelman and Toda in motivic
Donaldson–Thomas theory.

###### Contents

  1. 1 Introduction
    1. 1.1 Main result

##  1\. Introduction

###  1.1. Main result

Throughout this article, we will be working over the field of complex numbers
$\basefield$. Motivated by applications to Donaldson–Thomas invariants (see
LABEL:motivations), this paper studies the gluing of local invariants attached
to a $(-1)$–shifted symplectic derived Deligne–Mumford stack $X$ [1111.3209].
The simplest example of such $X$ is the derived critical locus $\dCrit(U,f)$
of a function $f$ on a smooth scheme $U$ (see
LABEL:symmetryderivedcriticalocusexplicit) and a Darboux lemma proved in
[MR3904157] shows that in fact, Zariski locally, any $(-1)$–shifted symplectic
$X$ is symplectically equivalent to some $\dCrit(U,f)$. This paper starts with
the observation that such local presentations are not unique:

######  Example 1.1.

Given the two pairs $(\affineline{\basefield},x^{2})$ and
$(\affine^{2}_{\basefield},x^{2}+y^{2})$, we have a canonical symplectic
identification
$\dCrit(\affineline{\basefield},x^{2})=\dCrit(\affine^{2}_{\basefield},x^{2}+y^{2})=\Spec(\basefield)$,
where $\Spec(\basefield)$ is equipped with the zero $(-1)$–shifted symplectic
form. More generally, we can always reproduce this phenomena by adding a
quadratic variable. For instance,
$\dCrit(\affineline{\basefield},x^{3})=\dCrit(\affine^{2}_{\basefield},x^{3}+y^{2})$.

Our goal is to investigate the moduli of such local Darboux parametrizations
of a $(-1)$–shifted symplectic derived Deligne–Mumford stack $X$ in a way that
resolves all ambiguities.

To start with, part of the moduli data should concern an _exact structure_ for
the symplectic form. Indeed, if $X=\dCrit(U,f)$ then the symplectic structure
has a canonical exact structure (see
LABEL:symmetryderivedcriticalocusexplicit). More generally, a result of
Deligne implies that every $(-1)$–shifted symplectic derived stack admits an
exact structure – see [MR3285853, Cor. 5.3] and our refinement in
LABEL:canonicalexactstructure. For this reason, we will consider
$(-1)$–shifted symplectic derived Deligne–Mumford stack $X$ _together with a
choice of exact structure_.

The naive attempt to define a moduli functor would be to assign to every
affine Zariski open $S=\Spec(A)$ in $X$, the $\infty$–groupoid classifying
_exact symplectic identifications_ of $S$111equipped with the exact symplectic
structure restricted from $X$. with a derived critical locus $\dCrit(U,f)$ for
$U$ a smooth scheme and $f$ a function on $U$. Unfortunately, this naive
attempt is not even functorial: if $S^{\prime}\subseteq S$ is an open Zariski
subset of $S$, there is a priori no canonical way to produce a Zariski open
$U^{\prime}$ of $U$ with a function $f^{\prime}$ whose derived critical locus
is $S^{\prime}$. In order to overcome this difficulty we observe that
$\dCrit(U,f)$ with its exact $(-1)$–symplectic structure only depends on the
smooth formal scheme given by the formal completion $\formalU$ of the derived
critical locus $\dCrit(U,f)$ inside $U$ together with the Taylor expansion of
$f$ (see LABEL:dCritonlydependsformal). Therefore, we can replace _algebraic_
smooth schemes $U$, by _formal_ smooth schemes $\formalU$ in our Darboux local
models. The fact that étale covers can be uniquely extended along
infinitesimal thickenings solves the functoriality problem.

Following this discussion, we can introduce the central object of this work,
namely, the Darboux moduli


# Matters arising:
Is there evidence of a trend in the CO2 airborne fraction?

Mikkel Bennedsen, Eric Hillebrand, Siem Jan Koopman  Department of Economics
and Business Economics and CREATES, Aarhus University, Fuglesangs Allé 4, 8210
Aarhus V, Denmark. E-mail<EMAIL_ADDRESS>Department of Economics and
Business Economics and CREATES, Aarhus University, Fuglesangs Allé 4, 8210
Aarhus V, Denmark. E-mail<EMAIL_ADDRESS>Department of Econometrics,
School of Business and Economics, Vrije Universiteit Amsterdam, De Boelelaan
1105, 1081 HV Amsterdam, The Netherlands, and CREATES. E-mail:
<EMAIL_ADDRESS>
In a paper recently published in this journal, van Marle et al. (van Marle et
al., 2022) introduce an interesting new data set for land use and land cover
change CO2 emissions (LULCC) that they use to study whether a trend is present
in the airborne fraction (AF), defined as the fraction of CO2 emissions
remaining in the atmosphere. Testing the hypothesis of a trend in the AF has
attracted much attention, with the overall consensus that no statistical
evidence is found for a trend in the data (Knorr, 2009; Gloor et al., 2010;
Raupach et al., 2014; Bennedsen et al., 2019). In their paper, van Marle et
al. analyze the AF as implied by three different LULCC emissions time series
(GCP, H&N, and their new data series). In a Monte Carlo simulation study based
on their new LULCC emissions data, van Marle et al. find evidence of a
declining trend in the AF.

In this note, we argue that the statistical analysis presented in van Marle et
al. can be improved in several respects. Specifically, the Monte Carlo study
presented in van Marle et al. is not conducive to determine whether there is a
trend in the AF. Further, we re-examine the evidence for a trend in the AF by
using a variety of different statistical tests. The statistical evidence for
an uninterrupted (positive or negative) trend in the airborne fraction remains
mixed at best. When allowing for a break in the trend, there is some evidence
for upward trends in both subsamples.

## The Monte Carlo design of van Marle et al.

Monte Carlo studies can be employed in situations where a test-statistic has
unknown properties (size and/or power), either because it is suspected or
known that the data-generating process (DGP) does not satisfy the assumptions
of the test, or because a relatively small sample size makes the invocation of
a central limit theorem for the test-statistic dubious. In a Monte Carlo
design, the DGP needs to be well enough understood to specify a realistic and
meaningful set of processes to simulate from. To establish the size of the
test, $Pr(\textrm{reject }H_{0}|H_{0}\textrm{ is true})$, the DGP is simulated
under the null hypothesis (no trend), and the test-statistic is calculated for
a large number of simulated trajectories. The percentiles of the ensemble of
simulated test-statistics reveal how one should choose a critical value for
the test in order to obtain a desired size. One can compare the test-statistic
obtained from the data to these percentiles in order to estimate the
probability of observing the result from the data if the null hypothesis is
true. To determine the power of a test, $Pr(\textrm{reject
}H_{0}|H_{1}\textrm{ is true})$, the DGP is simulated under the alternative
hypothesis (trend is present). This necessitates a defendable choice of the
trend specification for the simulated DGP. The percentiles of the ensemble of
simulated test-statistics reveal how one should choose a critical value for
the test in order to obtain a desired power.

In their paper, Marle et al. perturb the original data with simulated Gaussian
random variables, where variances are chosen that correspond to empirical
variances of the data time series. The non-parametric Mann-Kendall statistic
is then computed to test for a trend in the simulated AF time series. These
two steps are repeated $10,000$ times, and the number of times a positive or
negative trend is detected by the Mann-Kendall test on the simulated data is
counted (Table 1 in van Marle et al.). While this design has some elements of
the procedure outlined above, it falls short of giving the researcher a handle
on whether they are simulating under the null hypothesis of no trend or under
the alternative hypothesis. Therefore, the frequency of detected trends in the
design of van Marle et al. cannot be interpreted as a probability under either
hypothesis. Similarly, it cannot be interpreted as a probability that the
original data contain a trend, as appears to be suggested in van Marle et al.

## How much statistical evidence is there for a negative trend?

There are many different ways to test for the presence of a trend, just as
there are many conceptualizations of trends (Hamilton, 1994). In the Monte
Carlo simulation study of van Marle et al., the non-parametric Mann-Kendall
statistic is used to test for a trend. However, the Mann-Kendall test-
statistic is not computed for the actual data set. When applied to the data,
no statistical evidence of a trend is found (Table 1).

Table 1: $p$-values for various tests on the AF data from van Marle et al. The $p$-values are reported for the Mann-Kendall tests ( $H_{0}^{MK}:\textnormal{``data do not contain a trend''}$ against the two-sided alternative $H_{1}:\textnormal{``trend in data''}$; $H_{1}^{MK+}$: the one-sided alternative for “positive trend in data”; $H_{1}^{MK-}$: the one-sided alternative for “negative trend in data”), for the trend slope coefficients ( $H_{0}^{\textrm{slope1}}:b=0$ against the two-sided alternative $H_{1}:b\neq 0$ for the linear trend model, $AF_{t}=a+bt+\epsilon_{t}$; $H_{0}^{\textrm{slope2}}:b=0$ against the two-sided alternative $H_{1}:b\neq 0$ for the linear trend model with intercept break in 1988 for raw series and in 1990 for filtered series, $AF_{t}=a_{1}+bt+a_{2}I(t\geq\tau)+\epsilon_{t}$) and for the break coefficients ( $H_{0}^{\textrm{br-trend}}:b_{2}=0$ against the two-sided alternative $H_{1}:b_{2}\neq 0$ for the linear trend model with break, $AF_{t}=a_{1}+b_{1}t+a_{2}I(t\geq\tau)+b_{2}(t-\tau)I(t\geq\tau)+\epsilon_{t}$; $H_{0}^{\textrm{br-intercept}}:a_{2}=0$ against the two-sided alternative $H_{1}:a_{2}\neq 0$ for the linear trend model with break, $AF_{t}=a_{1}+bt+a_{2}I(t\geq\tau)+\epsilon_{t}$), for the three different data series studied in van Marle et al. (GCP, H$\&$N, New), using both the raw and filtered versions (‘raw’ and ‘filter’). |  | Mann-Kendall tests |  | Regression tests
---|---|---|---|---
|  | $H_{0}^{MK}$ | $H_{0}^{MK+}$ | $H_{0}^{MK-}$ |  | $H_{0}^{\textrm{slope1}}$ | $H_{0}^{\textrm{slope2}}$ | $H_{0}^{\textrm{br-trend}}$ | $H_{0}^{\textrm{br-intercept}}$
GCP | raw | 0.1017 | 0.0508 | 0.9492 |  | 0.1447 | 0.0174 | 0.1791 | 0.0456
| filter | 0.0646 | 0.0323 | 0.9677 |  | 0.0634 | 0.0030 | 0.2254 | 0.0213
H&N | raw | 0.4819 | 0.2409 | 0.7591 |  | 0.8262 | 0.0184 | 0.3514 | 0.0096
| filter | 0.6230 | 0.3115 | 0.6885 |  | 0.8738 | 0.0048 | 0.7698 | 0.0017
New | raw | 0.7415 | 0.6292 | 0.3708 |  | 0.3473 | 0.0184 | 0.3692 | 0.0096
| filter | 0.1850 | 0.9075 | 0.0925 |  | 0.1130 | 0.0474 | 0.8109 | 0.0016

To complement the Mann-Kendall test, we perform a number of further tests for
a trend in the AF, none of them finding evidence of a declining trend in the
AF (Supplementary Information). Here, we report the results for the test that
most researchers would resort to at first pass: the $t$ test for for the
least-squares estimate of a linear trend in the linear model
$AF_{t}=a+bt+\epsilon_{t}$. The estimated slope coefficient $\hat{b}$ is
positive for the GCP and H&N data series, while for the AF series computed
from the new LULCC series, the slope estimate is indeed negative
(Supplementary Information, Table 2). These statements hold for both the raw
and the filtered series. Figure 1 presents the fitted trends (blue dashed
line). The corresponding tests indicate that the slope estimates are not
significant for any of the AF time series ($H_{0}^{\textrm{slope1}}$, Table
1).

Breaks in the level and trend of the AF have been discussed in Keenan et al.
(Keenan et al., 2016). The time series econometrics literature has a tradition
of studying these “structural breaks” (Perron, 2006; Clements and Hendry,
2006; Mills, 2010). We consider the single-break model (Perron and Zhu, 2005):

$AF_{t}=a_{1}+b_{1}t+a_{2}I(t\geq\tau)+b_{2}(t-\tau)I(t\geq\tau)+\epsilon_{t},$

where $\tau\geq 1$ and $\epsilon_{t}$ is a covariance-stationary error. This
model allows for a change in the intercept from $a_{1}$ to $a_{1}+a_{2}$ and a
change in the trend slope from $b_{1}$ to $b_{1}+b_{2}$ at time $\tau$. The
coefficients are estimated for each $\tau=11,\ldots,51$ (leaving out 10
observations at each side of the sample to avoid degeneracies). The break
point $\tau$ is chosen where the sum of squared errors is lowest.

This method finds break points in the three raw series in 1988 and in the
three filtered series in 1990. It is not without peril to connect
statistically found change points to historical events. In the case of the AF
data, relevant events may be the Pinatubo minimum in the change in atmospheric
concentrations (Bousquet et al., 2000; Angert et al., 2004), the collapse of
the Soviet Union and the first Gulf war in anthropogenic emissions, and fire
events documented for this period in the LULCC series (van Marle et al.,
2022). We can subsequently test for whether the data support a break in the
trend, by testing the hypothesis $b_{2}=0$. For all data series, this
hypothesis cannot be rejected at a nominal level of $5\%$ ($H_{0}^{\textrm{br-
trend}}$, Table 1), indicating that there is no break in the trend.
Conversely, the hypothesis $a_{2}=0$ can be rejected at a nominal level of
$5\%$ ($H_{0}^{\textrm{br-intercept}}$, Table 1), indicating that the break is
driven by a change in the intercept. Figure 1 presents the resulting estimated
trends (red solid line). All segmented trend estimates are positive
(Supplementary Information, Table 5) and highly significant
($H_{0}^{\textrm{slope2}}$, Table 1).

Figure 1: AF data (black, solid), estimated trends from the linear trend model
(blue, dashed), and the linear trend model with a break (red, solid).

## References

  * Angert et al. (2004) Angert, A., S. Biraud, C. Bonfils, W. Buermann, and I. Fung (2004).  CO2 seasonality indicates origins of post-Pinatubo sink.  Geophysical Research Letters 31(11), L11103.
  * Bennedsen et al. (2019) Bennedsen, M., E. Hillebrand, and S. J. Koopman (2019).  Trend analysis of the airborne fraction and sink rate of anthropogenically released CO2.  Biogeosciences 16(18), 3651–3663.
  * Bousquet et al. (2000) Bousquet, P., P. Peylin, P. Ciais, C. Le Quéré, P. Friedlingstein, and P. P. Tans (2000).  Regional changes in carbon dioxide fluxes of land and oceans since 1980\.  Science 290(5495), 1342–1346.
  * Clements and Hendry (2006) Clements, M. P. and D. F. Hendry (2006).  Forecasting with breaks.  In G. Elliott, C. W. J. Granger, and A. Timmermann (Eds.), Handbook of Economic Forecasting: Volume 1, pp. 605–657. Elsevier.
  * Durbin and Koopman (2012) Durbin, J. and S. J. Koopman (2012).  Time series analysis by state space methods.  Number 38. Oxford University Press.
  * Fuller (1996) Fuller, W. A. (1996).  Introduction to Statistical Time Series.  John Wiley & Sons.
  * Gloor et al. (2010) Gloor, M., J. L. Sarmienti, and N. Gruber (2010).  What can be learned about carbon cycle climate feedbacks from the CO2 airborne fraction?  Atmospheric Chemistry and Physics 10, 7739 – 7751.
  * Hamilton (1994) Hamilton, J. D. (1994).  Time Series Analysis.  Princeton University Press.
  * Harvey (1989) Harvey, L. D. D. (1989).  Managing atmospheric CO2.  Climatic Change 15(3), 343–381.
  * Johansen (1995) Johansen, S. (1995).  Likelihood-Based Inference in Cointegrated Vector Autoregressive Models.  Oxford University Press.
  * Keenan et al. (2016) Keenan, T. F., I. C. Prentice, J. G. Canadell, C. A. Williams, H. Wang, M. Raupach, and G. J. Collatz (2016).  Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake.  Nature communications 7(1), 1–10.
  * Knorr (2009) Knorr, W. (2009).  Is the airborne fraction of anthropogenic CO2 emissions increasing?  Geophysical Research Letters 36.
  * Mills (2010) Mills, T. C. (2010).  “Skinning a cat”: alternative models of representing temperature trends.  Climatic Change 101(3), 415–426.
  * Newey and West (1987) Newey, W. K. and K. D. West (1987).  A simple, positive semi-definite, heteroskedasticity and autocorrelationconsistent covariance matrix.  Econometrica 55(394), 703–708.
  * Perron (2006) Perron, P. (2006).  Dealing with structural breaks.  In K. Patterson and T. C. Mills (Eds.), Palgrave Handbook of Econometrics, Vol. 1: Econometric Theory, pp. 278–352. Palgrave Macmillan.
  * Perron and Zhu (2005) Perron, P. and X. Zhu (2005).  Structural breaks with deterministic and stochastic trends.  Journal of Econometrics 129(1-2), 65–119.
  * Raupach et al. (2014) Raupach, M. R., M. Gloor, J. L. Sarmiento, J. G. Canadell, T. L. Frölicher, T. Gasser, R. A. Houghton, C. Le Quéré, and C. M. Trudinger (2014).  The declining uptake rate of atmospheric CO2 by land and ocean sinks.  Biogeosciences 11(13), 3453–3475.
  * van Marle et al. (2022) van Marle, M. J. E., D. van Wees, R. A. Houghton, R. D. Field, J. Verbesselt, and G. R. van der Werf (2022).  New land-use-change emissions indicate a declining CO2 airborne fraction.  Nature 603(7901), 450–454.

##  Appendix A Supplementary Information

###  A.1 Evidence from linear trend regressions, with and without breaks

We consider the linear trend regression model

$\displaystyle AF_{t}=a+bt+\epsilon_{t},$ (A.1)

where $\epsilon_{t}$ is an error term and $AF_{t}$ is one of the measures of
AF (e.g. the ‘raw’ or ‘filtered’ AF). The purpose is to statistically test the
null hypothesis of “no trend”, formally given as

$\displaystyle H_{0}:b=0,$

against one of the following alternative hypotheses

$\displaystyle H_{0}^{two-sided}:b\neq 0,$ $\displaystyle H_{0}^{pos}:b>0,$
$\displaystyle H_{0}^{neg}:b<0.$

We can perform tests on the parameter $b$ using the least squares method. We
consider the $t$ test-statistic

$\displaystyle t=\frac{\hat{b}}{s.e(\hat{b})},$

where $\hat{b}$ is the least squares estimator of $b$ from Equation (A.1) and
$s.e(\hat{b})$ is the standard error of this estimator, which we have
calculated using the heteroskedasticity and autocorrelation robust (HAC)
estimator of $s.e.(\hat{b})$ proposed by Newey and West (1987). We also
present the results obtained from using the standard estimator of
$s.e(\hat{b})$ which lead to similar findings and conclusions.

Table 2 presents the regression estimation results obtained from the least
squares method applied to Equation (A.1), while Table 3 presents the
$p$-values from the test of $H_{0}:b=0$. From the results in Table 3 we learn
that it is not possible to conclude that there is a trend in the data (either
positive or negative), no matter which data set is used.

Table 2: Least-squares regression output of the linear model. Output of the least-squares regression of the model $y_{t}=a+bt+\epsilon_{t}$, where $t=0,1,\ldots,60$. Standard errors are calculated using the Newey and West (1987) HAC estimator. Data: | GCP-raw | GCP-filter | H&N-raw | H&N-filter | New-raw | New-filter
---|---|---|---|---|---|---
$\hat{a}$ (intercept) | 0.4108 | 0.4113 | 0.4657 | 0.4662 | 0.5120 | 0.5132
SE($\hat{a}$) | 0.0308 | 0.0213 | 0.0340 | 0.0211 | 0.0363 | 0.0226
$t$-stat | 13.3436 | 19.3437 | 13.7112 | 22.1169 | 14.1083 | 22.6800
$\hat{b}$ (trend) | 0.0011 | 0.0010 | 0.0002 | 0.0001 | -0.0008 | -0.0009
SE($\hat{b}$) | 0.0007 | 0.0005 | 0.0008 | 0.0005 | 0.0009 | 0.0006
$t$-stat | 1.4585 | 1.8567 | 0.2195 | 0.1589 | -0.9398 | -1.5847
Table 3: Testing for a trend in AF using least squares and HAC standard errors. Testing $H_{0}:b=0$ using least squares and HAC standard errors. We have carried out the test for the three different data series studied in the main paper (GCP, H&N, New), as well as using the raw and filtered data (‘raw’ and ‘filter’). |  $p$-value ($H_{0}^{two-sided}$) |  $p$-value ($H_{0}^{pos}$) |  $p$-value ($H_{0}^{neg}$)
---|---|---|---
GCP (raw) | 0.1447 | 0.0724 | 0.9276
GCP (filter) | 0.0634 | 0.0317 | 0.9683
H&N (raw) | 0.8262 | 0.4131 | 0.5869
H&N (filter) | 0.8738 | 0.4369 | 0.5631
New (raw) | 0.3473 | 0.8263 | 0.1737
New (filter) | 0.1130 | 0.9435 | 0.0565

Next, the focus is on the statistical evidence of a trend in AF, after
allowing for a break in the level and in the slope of AF. We consider the
single-break model (Perron and Zhu, 2005):

$\displaystyle
AF_{t}=a_{1}+a_{2}I(t\geq\tau)+b_{1}t+b_{2}tI(t\geq\tau)+\epsilon_{t},$ (A.2)

where $\epsilon_{t}$ is a covariance-stationary error. This model allows for a
change in the intercept from $a_{1}$ to $a_{1}+a_{2}$ and a change in the
trend slope from $b_{1}$ to $b_{1}+b_{2}$ at time $\tau$. The coefficients are
estimated for each $\tau=11,\ldots,51$ (leaving out 10 observations at each
side of the sample to avoid degeneracies). The break point $\tau$ is chosen
where the sum of squared errors is lowest. This method detects single-break
points in the three raw series in 1988 and in the three filtered series in
1990.

Table 4 presents the estimation results obtained from the least squares method
applied to Equation (A.2); we find that the estimates of $b_{2}$ are not
significantly different from zero. Hence, we conclude that the data prefer a
model with no break in the slope of the trend.

Table 4: Least-squares regression output of the linear model with break in both intercept and trend. Output of the least-squares regression of the model $AF_{t}=a_{1}+b_{1}t+a_{2}\cdot I(t\geq\tau)+b_{2}\cdot(t-\tau)\cdot I(t\geq\tau)+\epsilon_{t}$, where $t=0,1,\ldots,60$ and $\tau=30$ for the raw data series (corresponding to a break in year $1988$) and $\tau=32$ for the filtered data series (corresponding to a break in year $1990$). Standard errors are calculated using the Newey and West (1987) HAC estimator. Data: | GCP-raw | GCP-filter | H&N-raw | H&N-filter | New-raw | New-filter
---|---|---|---|---|---|---
$\hat{a}_{1}$ (intercept) | 0.3477 | 0.3755 | 0.4029 | 0.4327 | 0.4478 | 0.4785
SE($\hat{a}_{1}$) | 0.0390 | 0.0228 | 0.0411 | 0.0248 | 0.0425 | 0.0263
$t$-stat | 8.9245 | 16.4775 | 9.7948 | 17.4312 | 10.5381 | 18.1986
$\hat{b}_{1}$ (trend) | 0.0069 | 0.0040 | 0.0064 | 0.0033 | 0.0055 | 0.0024
SE($\hat{b}_{1}$) | 0.0022 | 0.0014 | 0.0025 | 0.0016 | 0.0026 | 0.0017
$t$-stat | 3.1870 | 2.9582 | 2.5457 | 2.0827 | 2.1010 | 1.3977
$\hat{a}_{2}$ (break, intercept) | -0.1462 | -0.0791 | -0.1856 | -0.1185 | -0.1921 | -0.1262
SE($\hat{a}_{2}$) | 0.0516 | 0.0288 | 0.0533 | 0.0299 | 0.0536 | 0.0313
$t$-stat | -2.8324 | -2.7444 | -3.4826 | -3.9649 | -3.5827 | -4.0324
$\hat{b}_{2}$ (break, trend) | -0.0042 | -0.0023 | -0.0031 | -0.0006 | -0.0031 | -0.0005
SE($\hat{b}_{2}$) | 0.0032 | 0.0019 | 0.0033 | 0.0021 | 0.0034 | 0.0022
$t$-stat | -1.3434 | -1.2124 | -0.9318 | -0.2927 | -0.8980 | -0.2392

These findings provide some support in considering the model

$\displaystyle y_{t}=a_{1}+bt+a_{2}\cdot I(t\geq\tau)+\epsilon_{t},$ (A.3)

where $t=0,1,\ldots,60$ and $\tau=30$ for the ‘raw’ data series (corresponding
to a break in year $1988$) and $\tau=32$ for the ‘raw’ data series
(corresponding to a break in year $1990$). Table 5 presents the estimation
results obtained from the least squares method applied to Equation (A.3).

Table 5: Least-squares regression output of the linear model with break in intercept. Output of the least-squares regression of the model $AF_{t}=a_{1}+bt+a_{2}\cdot I(t\geq\tau)+\epsilon_{t}$, where $t=0,1,\ldots,60$ and $\tau=30$ for the ‘raw’ data series (corresponding to a break in year $1988$) and $\tau=32$ for the ‘raw’ data series (corresponding to a break in year $1990$). Standard errors are calculated using the Newey and West (1987) HAC estimator. Data: | GCP-raw | GCP-filter | H&N-raw | H&N-filter | New-raw | New-filter
---|---|---|---|---|---|---
$\hat{a}$ (intercept) | 0.3817 | 0.3918 | 0.4278 | 0.4370 | 0.4726 | 0.4822
SE($\hat{a}_{1}$) | 0.0298 | 0.0184 | 0.0312 | 0.0198 | 0.0322 | 0.0210
$t$-stat | 12.8039 | 21.2541 | 13.7012 | 22.0220 | 14.6593 | 22.9481
$\hat{b}$ (trend) | 0.0045 | 0.0029 | 0.0046 | 0.0030 | 0.0038 | 0.0022
SE($\hat{b}$) | 0.0014 | 0.0009 | 0.0014 | 0.0010 | 0.0015 | 0.0011
$t$-stat | 3.2914 | 3.0982 | 3.1824 | 3.0596 | 2.5415 | 2.0714
$\hat{a}_{2}$ (break, intercept) | -0.1378 | -0.0791 | -0.1794 | -0.1185 | -0.1860 | -0.1262
SE($\hat{a}_{2}$) | 0.0519 | 0.0301 | 0.0529 | 0.0298 | 0.0532 | 0.0312
$t$-stat | -2.6531 | -2.6280 | -3.3905 | -3.9760 | -3.4993 | -4.0474

For completeness, we also compare the two models in (A.2) and (A.3) using
statistical goodness-of-fit tests. We compare the likelihoods of the two
models, evaluated at their optimized parameters and assuming that the error
term $\epsilon_{t}$ is independent, identically and normally distributed, by
means of the Bayesian Information Criteria (BIC). Since the model in (A.2)
nests the model in (A.3), it is necessarily the case that the former has a
higher likelihood than the latter. When looking at the BIC, however, we see
that the smaller model (A.3) has a lower BIC value than the larger model
(A.2), implying that the data prefer the smaller model. Finally, since the two
models are nested, we can also consider a likelihood ratio test and test the
hypothesis $H_{0}:b_{2}=0$ in Equation (A.2). Given the results presented in
Table 6, this null hypothesis cannot be rejected. These results confirm the
finding above: the data prefer a model with a break in the intercept, but no
break in the trend.

Table 6: Statistical comparison of the two models. Statistical comparison of the model with break in intercept (Model 1) and the model with break in intercept and trend (Model 2). Model 1: $AF_{t}=a_{1}+bt+a_{2}\cdot I(t\geq\tau)+\epsilon_{t}$. Model 2): $AF_{t}=a_{1}+b_{1}t+a_{2}\cdot I(t\geq\tau)+b_{2}\cdot(t-\tau)\cdot I(t\geq\tau)+\epsilon_{t}$. Here, $t=0,1,\ldots,60$ and $\tau=30$ for the ‘raw’ data series (corresponding to a break in year $1988$) and $\tau=32$ for the ‘filter’ data series (corresponding to a break in year $1990$). “logL” is the log-likelihood of the model obtained by least-squares regression, i.e. $logL=-\frac{1}{2}n\log(SSE/n)$, where $SSE$ is the sum of squared errors and $n=61$ is the number of observations. “BIC” is the Bayesian information criterion given by $BIC=-2logL+k\log n$, where $k$ is the number of parameters in a given model, so that $k=3$ in the case of Model 1 and $k=4$ in the case of Model 2. “LR-stat” is the likelihood ratio test-statistic for the null that the break in the trend is absent, i.e. $b_{2}=0$, that is $LR=-2(logL_{1}-logL_{2})$. Under the null, we have $LR\sim\chi^{2}(1)$, meaning that the critical value for the likelihood ratio (LR) test at a $5\%$ level is $3.8415$.  Data: | GCP-raw | GCP-filter | H&N-raw | H&N-filter | New-raw | New-filter
---|---|---|---|---|---|---
logL (Model 1) | 137.0926 | 168.0499 | 134.2751 | 163.5550 | 132.3149 | 160.0704
logL (Model 2) | 138.0433 | 168.8265 | 134.7362 | 163.6008 | 132.7434 | 160.1010
BIC (Model 1) | -261.8527 | -323.7672 | -256.2176 | -314.7773 | -252.2972 | -307.8082
BIC (Model 2) | -259.6432 | -321.2095 | -253.0289 | -310.7581 | -249.0433 | -303.7586
LR-stat | 1.9014 | 1.5531 | 0.9222 | 0.0916 | 0.8570 | 0.0612

###  A.2 Evidence from year-to-year changes

Year-to-year changes in a time series $y_{t}$ are measured as $\Delta
y_{t}=y_{t}-y_{t-1}$. The average change

$\frac{1}{T-1}\sum_{t=2}^{T}\Delta y_{t},$

is an indication of the direction of a trend, if there is one: a negative
average indicates a negative trend. We consider the running average change
defined as $(\tau-1)^{-1}\sum_{t=2}^{\tau}\Delta y_{t}$, for
$\tau=2,\ldots,T$, together with the 95% point-wise confidence intervals for
each change average (Figure 2). All presented averages indicate change
estimates which are not significantly different from zero.

Figure 2: Running average change with increasing sample of year-to-year
changes for the three different data series studied in the main paper (row-
wise: GCP, H&N, New), as well as using the raw and filtered data (column-wise:
‘raw’ and ‘filter’). Running average is defined as
$(\tau-1)^{-1}\sum_{t=2}^{\tau}\Delta y_{t}$, for $\tau=2,\ldots,T$.

###  A.3 Evidence from stochastic trends

Trends can be conceptualized as deterministic or as stochastic (Johansen,
1995; Fuller, 1996). A standard model for the simplest form of a stochastic
trend, a random walk, is the local level model (Harvey, 1989; Durbin and
Koopman, 2012). The local level model is also known as the random walk plus
noise model, and it is given by

$y_{t}=\mu_{t}+\varepsilon_{t},\qquad\mu_{t+1}=\mu_{t}+\eta_{t},\qquad
t=1,\ldots,T,$

where $\mu_{1}$ is treated as unknown. The estimates of the level $\mu_{t}$
(the stochastic trend) using all observations $y_{1},\ldots,y_{T}$ are
presented in Figure 3. The 95% point-wise confidence intervals are also
provided.

Figure 3: Level $\mu_{t}$ estimates, with 95% confidence interval, from the
local level model for the three different data series studied in the main
paper (row-wise: GCP, H&N, New), as well as using the raw and filtered data
(column-wise: ‘raw’ and ‘filter’).

The evidence for trends in all series is rather weak. The estimated stochastic
trends for the GCP and H&N data series (raw and filtered) show slowly upward
moving trends. The new AF series have estimated trends that are somewhat
downward moving. However, the patterns are not showing continuous declines.
These trends reveal that they are more or less constant from 1959 to about
1990, then shifts take place to lower levels, which remain more or less
constant, that is, from the mid-90s to the end of the sample. To follow-up on
these findings, we introduce a break dummy $I_{\\{t>\tau\\}}$ in the random
walk specification for the trend $\mu_{t}$, which takes value one for $t>\tau$
and zero elsewhere. We record the $t$ tests of the break dummy
$I_{\\{t>\tau\\}}$ for each $\tau=1,\ldots,T$ in Figure 4. There is little
evidence of a break in the GCP and H&N time series. There is, however, strong
evidence of a break in the AF series (raw and filtered) around the years
1988-1991. These findings confirm our earlier conclusions.

Figure 4: $t$ tests for break at time $\tau$ in local level estimates, with
$\tau=1,\ldots,T$, for the three different data series studied in the main
paper (row-wise: GCP, H&N, New), as well as using the raw and filtered data
(column-wise: ‘raw’ and ‘filter’).


# Holographic Abrikosov lattice: vortex matter from black hole

Chuan-Yin Xia1<EMAIL_ADDRESS>Hua-Bi Zeng1<EMAIL_ADDRESS>Yu
Tian2,3<EMAIL_ADDRESS>Chiang-Mei Chen4<EMAIL_ADDRESS>Jan Zaanen5
<EMAIL_ADDRESS>1Center for Gravitation and Cosmology, College of
Physical Science and Technology, Yangzhou University, Yangzhou 225009, China
2School of Physics, University of Chinese Academy of Sciences, Beijing 100049,
China  3Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing
100190, China  4Department of Physics, Center for High Energy and High Field
Physics (CHiP), National Central University, Chungli 32001, Taiwan  5Institute
Lorentz for Theoretical Physics, Leiden University, Leiden, The Netherlands

###### Abstract

The AdS/CFT correspondence provides a unique way to study the vortex matter
phases in superconductors. We solved the nonlinear equations of motion for the
Abelain-Higgs theory living on the AdS4 black hole boundary that is dual to a
two dimensional strongly coupled type II superconductor at temperature $T$
with a perpendicular external uniform magnetic field $B_{0}$. We found the
associated two critical magnetic fields, $B_{c1}(T)$ and $B_{c2}(T)$. For
$B_{0}<B_{c1}(T)$ the magnetic field will be expelled out by the
superconductor resembling the Meissner effect and the superconductivity will
be destroyed when $B_{0}>B_{c2}(T)$. The Abrikosov lattice appears in the
range $B_{c1}(T)<B_{0}<B_{c2}(T)$ including, due to the finite size and
boundary effect, several kinds of configurations such as hexagonal, square and
slightly irregular square lattices, when the magnetic field is increased. The
upper and lower critical fields behave as inverse squares of coherence length
and magnetic penetration depth respectively which matches the well known
consensus.

Introduction — A well-known property of the type II superconductors is the
quantization of the magnetic flux in the mixed state, where the magnetic field
penetrates the sample as vortices forming an Abrikosov lattice in which each
vortex carrying one thread of magnetic field with the flux $\Phi_{0}=hc/2e$
Abrikosov . Close to the transition temperature, the formation of vortex
lattice can be simulated by the time dependent Ginzburg-Landau (GL) equation
with parameter $\kappa>1/\sqrt{2}$. The GL parameter $\kappa$ is defined as
the ratio of the magnetic penetration depth $\lambda$ to the coherence length
$\xi$ of the order parameter, $\kappa=\lambda/\xi$, and it can be calculated
from the microscopic parameters of the material within the
Bardeen–Cooper–Schrieffer (BCS) theory. In a type I superconductor,
$\kappa<1/\sqrt{2}$, the interaction between vortices is purely attractive
which results in their fusion into macroscopic normal domains in the
intermediate state, then no stable vortices appear. On the contrary in a type
II superconductor, $\kappa>1/\sqrt{2}$, the interaction between vortices is
purely repulsive so the vortices are stable and form a lattice in the mixed
state Brandt ; Bogomolnyi1 ; Bogomolnyi2 ; Jacobs ; Beekman:2011tn . For
review of GL theory of type II superconductors under magnetic field, please
refer to Ginzburg ; Gennes ; Saint ; Blatter:1994zz ; Tinkham ; Rosenstein .

In the AdS/CFT correspondence framework Maldacena ; Gubser ; Witten , the
holographic version of the superconducting model was firstly proposed in
Gubser2008 ; Hartnoll . The holographic superconductor model includes a
charged scalar field living in an AdS planar black hole. The scalar field,
when the temperature of the black hole is lower enough, can have a nonzero
profile with lower free energy than the trivial zero solution. Such an $U(1)$
symmetry broken configuration is dual to a superconducting state in the
boundary field theory by computing the conductivity via the AdS/CFT
correspondence dictionary Hartnoll . The computation of the associated GL
parameter $\kappa$ indicates that such a superconductor is always of type II
in the probe limit Umeh:2009ea ; Hartnoll:2008kx ; huabi . Afterward many
efforts have been devoted to consider external magnetic field effects
Nakano:2008xc ; Albash:2008eh ; Albash:2009ix , in particular, to find the
stable vortex state in the $s$-wave holographic superconductor, by solving the
time independent equation of motion (EoMs) for the scalar and gauge fields in
the bulk spacetime Montull ; Domenech ; Keranen ; Albash:2009iq ; Maeda ;
Aristomenis . Due to the highly nonlinear properties of the EoMs, single
vortex solution, rather than vortex lattice, had been obtained Montull ;
Domenech ; Keranen ; Albash:2009iq , while a static vortex lattice solution
was obtained by the perturbative method Maeda ; Aristomenis .

Rather than solving the time independent EoMs, the vortex lattice states can
also be obtained by studying the dynamics of a homogeneous superconductor in
the presence of a uniform magnetic field by solving the full time dependent
EoMs of the holographic superconductor model. The equilibrium phase within
fixed temperature and magnetic field can be obtained by finding the final
stable configuration which does not change in time anymore. This is very
similar to the vortex lattice formation dynamics simulation by solving the
time dependent Gindzburg-Landau equation Dorsey ; Du1 . In this letter, we
report the results of the equilibrium vortex lattice phases, the magnetization
curve and two critical fields in the phase diagram. We also address
temperature dependence of the magnetic penetration depth, the coherence length
and two critical fields which match the experimental observations.

Holographic Model — The action of the Abelian-Higgs model in AdS4 black holes,
in the unit $\hbar=c=G_{N}=1$, reads Gubser2008 ; Hartnoll

$S=\int
d^{4}x\sqrt{-g}\left(-\frac{1}{4}F^{2}-|D\Psi|^{2}-m^{2}|\Psi|^{2}\right),$
(1)

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and
$D_{\mu}=\nabla_{\mu}-iqA_{\mu}$ with the charge of scalar field, as a Cooper
pair, $q=2e$. The AdS4 black hole background in the Eddington-Finkelstein
coordinates is

$ds^{2}=\frac{\ell^{2}}{z^{2}}\left(-f(z)dt^{2}-2dtdz+dx^{2}+dy^{2}\right),$
(2)

in which $\ell$ is the AdS radius, $z$ is the radial coordinate of the AdS
bulk and $f(z)=1-(z/z_{h})^{3}$. Thus, $z=0$ is the AdS boundary while
$z=z_{h}$ is the horizon. The dual field theory lives at $z=0$, and the
information needed for the dual superconductor can be read from the behaviors
of the fields on the boundary by solving the dynamic coupled equations of
motion for $\Psi$ and $A_{\mu}$

$(D^{2}-m^{2})\Psi=0,\quad\nabla_{\mu}F^{\mu\nu}=i\Psi^{*}D^{\nu}\Psi-i\Psi(D^{\nu}\Psi)^{*}.$
(3)

Implicitly, the spontaneous broken of the local $U(1)$ symmetry for the field
theory is induced by a nonzero expectation value $\Psi^{(2)}$ of the scalar
operator dual to $\Psi$ in the bulk, which is read from the asymptotic
behavior of $\Psi$ on the boundary

$\Psi(z\sim 0,t,x,y)\approx\Psi^{(1)}(t,x,y)\,z+\Psi^{(2)}(t,x,y)\,z^{2},$ (4)

where the source $\Psi^{(1)}$ is set to be zero as a boundary condition when
solving the model. Furthermore, in order to introduce a magnetic field in the
dual holographic superconductor, the gauge fields on the boundary should be
dynamic. With the gauge fixing $A_{z}=0$, the behavior of the gauge fields on
the boundary is

$A_{\mu}(z\sim 0,t,x,y)\approx a_{\mu}(t,x,y)+b_{\mu}(t,x,y)\,z,$ (5)

in which $a_{\mu}$ can be regarded as the gauge field of the boundary theory,
while $b_{\mu}$ is related to the current $j_{\mu}$ as
$j_{\mu}=-b_{\mu}-\partial_{\mu}a_{t}+\partial_{t}a_{\mu}$ in the Eddington
coordinate following the holographic dictionary. We control the charge density
$\rho=-b_{t}$ in equivalent to turn the temperature. In the superconductor
case, we fix $j_{x}=j_{y}=0$ as the Neumann boundary condition for $A_{x}$ and
$A_{y}$ at $z=0$. Instead, in the superfluid case, the Dirichlet boundary
condition $a_{x}=a_{y}=0$ was imposed and the vortex lattice solution has been
found in Dias ; zeng ; XL ; Wittmer ; Ewerz ; Ankur .

Similar to the experimental setup for generating vortices, we prepare a
homogeneous superconducting state as the initial configuration, and then an
uniform external magnetic field is applied suddenly to the sample at $t=0$ by
turning on $A_{x}(t=0,z,x,y)=-B_{0}y/2$ and $A_{y}(t=0,z,x,y)=B_{0}x/2$. We
firstly prepare the initial homogeneous superconducting state at a fixed
temperature by the Newton–Raphson method. Its evolution under an external
magnetic field can be simulated by combining a Runge-Kutta method in the time
direction and a Chebyshev spectral method for the other three coordinates
$z,x,y$, similar to previous work on vortex lattice formation in a rotating
holographic superfluid zeng .

Figure 1: A hexagonal lattice at $T=0.95\,T_{c},\;B_{0}=0.428$: configurations
of (a) the order parameter, (b) the magnetic field, and (c) the radial
profiles of the order parameter and magnetic field in a single vortex.

Abrikosov Lattices — For the numerical simulation, we chose $\ell=1,m^{2}=-2$
and $q=1$ (implying $\Phi_{0}=2\pi$). In Fig. 1 we present a typical hexagonal
Abrikosov lattice solution as a final stable state does not change any more
for a sufficient long time in the dynamic process, when the temperature is
very close to $T_{c}$, and the configuration of order parameter and magnetic
field for one single vortex. Widths of the flux lines $\lambda$ and of the
order parameter defects $\xi$ can be fitted from the profile of magnetic field
$B(r)\sim 0.3949\exp(-r/\lambda)$ and the expectation value of the order
parameter $\langle O(r)\rangle\sim 1.3693\tanh(r/\sqrt{2}\xi)$. From Fig. 1,
we can estimate the values $\lambda\sim 1.579$ and $\xi\sim 1.1$,
respectively. Thus, the GL parameter is $\kappa\sim 1.435$, which belongs to
type II superconductors.

According to the GL theory analysis, the lattices with equilateral triangles
admit a slightly lower free energy than the square ones. It is interesting
that this result agrees with that of a simple argument based on the fact that
the triangular array is a “closed-packed” one, in which each vortex is
surrounded by a hexagonal array of other vortices. In this array, the nearest
neighbor distance can be evaluated from the averaged value of the magnetic
field in a vortex $\langle B\rangle$ as

$a_{\triangle}=\left(\frac{4}{3}\right)^{\frac{1}{4}}a_{\square}\approx\left(\frac{4}{3}\right)^{\frac{1}{4}}\left(\frac{\Phi_{0}}{\langle
B\rangle}\right)^{\frac{1}{2}}.$ (6)

Thus, for a given flux density, $a_{\square}<a_{\triangle}$. Taking into
account the mutual repulsion of the vortices, it is reasonable that the
structure with the greatest separation of the nearest neighbors would be
favored. From Eq. (6) the distance between two nearest vortices can be
computed as $a_{\triangle}\approx 1.075\sqrt{2\pi/0.178}\approx 6.39$, close
to the numerical simulation $a_{\triangle}\approx 6.63$.

Figure 2: Two representative vortex lattice solutions: (a-b)
$T=0.9\,T_{c},\;B_{0}=0.8$, and (c-d) $T=0.82\,T_{c},\;B_{0}=1.7$.

However, different things happen when the vortex number is increased by a
larger magnetic field. Two typical examples are given in Fig. 2. For the case
with $T=0.9\,T_{c},B_{0}=0.8$, there are 9 vortices forming a square pattern.
In this case, the distance between vortices, $a_{\square}\approx 5.02$, close
to the value from (6) $a_{\square}\approx 5.05$, is not small enough and the
finite size of system prohibits the formation of a hexagon pattern. For the
case $T=0.82\,T_{c},B_{0}=1.7$, the vortex number is 30, the distance between
two nearest vortices computed by (6) $a_{\triangle}\approx 3.43$ is also close
to the numerical simulation $a_{\triangle}\approx 3.214$. The finite size
effect is moderated, thus the hexagonal pattern is favored. However, the array
does not admit a perfect lattice configuration, which should be a consequence
of boundary effects since the vortices are close to the boundary. Keep
increasing the magnetic field to the $B_{c2}$, many superconducting areas
undergo a phase transition to normal state, leaving a superconducting island
with vortices crowded together without a well defined nearest vortex distance.

One can also clearly observe that the vortex size enlarges when the
temperature is increasing. More precisely, from the size of the vortex, and
similarly of the magnetic field profile, we can read out the dependence of
$\lambda$ and $\xi$ with respect to $T/T_{c}$ as shown in Fig. 3. Their
behaviors are, almost independent on $B_{0}$, consistent with the results by
GL theory Tinkham  for $T\approx T_{c}$

$\xi\sim
0.74\xi_{0}(1-T/T_{c})^{-1/2},\quad\lambda\sim\frac{\lambda_{0}}{\sqrt{2}}(1-T/T_{c})^{-1/2},$
(7)

with $\xi_{0}\sim 0.3313$ and $\lambda_{0}\sim 0.4825$. Thus the GL parameter
generally is $\kappa\sim 1.3916$. It is worth to note that these formulas can
fit the data for a broad range of temperature away from $T_{c}$, for example
with about 3% variation to the value at the point $T=0.95T_{c}$.

Figure 3: The temperature dependence of $\xi,\lambda$ with fitting formulae in
(7).

Magnetisation — To distinguish the Meissner phase and the vortex lattice
phase, we can compute the magnetisation $M$, which is defined as

$M(B_{0})=B_{0}-\left\langle B(x,y)\right\rangle,$ (8)

where the $B_{0}$ is the value of the external magnetic field applied at the
initial time, $B(x,y)$ is the magnetic field distribution in the final
equilibrium state. When the added external field increases from zero to
$B_{c1}$, the magnetic field is completely excluded then there is no magnetic
field inside the sample, i.e. $M=B_{0}$. While above $B_{c2}$, the
superconductivity is already completely destroyed, therefore $B(x,y)=B_{0}$
and $M$ should be zero. In the mixed state $B_{c1}<B_{0}<B_{c2}$, the
magnetization $M$ will decrease from $B_{c1}$ to zero gradually. In Fig. 4 we
show the magnetization versus $B_{0}$ for three different temperatures.

Figure 4: The magnetisation $M$ defined in Eq. (8) for three different
temperatures $T=0.82\,T_{c},0.9\,T_{c}$ and $0.95\,T_{c}$, from the curves we
can read the two critical fields where the magnetisation reaches its maximal
value and reduces to zero respectively.

Phase Diagram — From the magnetization curves we are able to read the two
critical magnetic fields $B_{c1}$ and $B_{c2}$ then to obtain the phase
diagram which is shown in Fig. 5. The upper critical field can be naively
estimated. As the magnetic field increases, more vortices enter and the
lattice becomes more compressed. At a certain point, the cores of vortices
overlap and no superconducting path is left for a transport current. Indeed,
the second critical field in the holographic superconductor model is

$B_{c2}\approx 16.64\,(1-T/T_{c}),$ (9)

which, according to the result (7), confirms the relation derived by the GL
theory Blatter:1994zz ; Tinkham

$B_{c2}=\frac{\Phi_{0}}{2\pi\xi^{2}}=\xi^{-2},$ (10)

where $\pi\xi^{2}$ is the size of the Abrikosov unit cell. Moreover, the lower
critical magnetic field $B_{c1}$ is close to the intuitive estimation for the
moment when the first vortex was created

$B_{c1}\approx\frac{\Phi_{0}}{2\pi\lambda^{2}}=\lambda^{-2},$ (11)

where $\pi\lambda^{2}$ is the magnetic field penetrated area.

Figure 5: The phase diagram of holographic superconductor under external
magnetic field.

Summary — Before the advent of AdS/CFT correspondence and the holographic
superconductor model defined in the AdS/CFT correspondence framework, the
study of the Abrikosov lattice formation dynamics was mainly based on the
Gindzburg-Landau theory. We find that the holographic superconductor model
offers another approach that the dynamics of magnetic quantum fluxes in
spatial two dimensions can be captured by solving the highly nonlinear coupled
PDEs in the bulk geometry. All the results agree with GL theory and the
experimental observations. There are many other issues needed to be studied,
for example, a detailed study of vortex matter near $B_{c2}$ may find the
vortex liquid state, extending the model to AdS5 will enable us to study the
vortex lines dynamics in high temperature superconductors. Also we focus on
the final equilibrium vortex lattice configuration in the present work, the
study of detailed vortex formation dynamics is in progress.

Acknowledgements — This work is supported by the National Natural Science
Foundation of China under Grant No. 11675140 (CYX, HBZ), 11975235, 12035016
(YT), and the Ministry of Science and Technology of the R.O.C. under the
grants MOST 109-2112-M-008-010, 110-2112-M-008-009 (CMC).

## References

  * (1) A. A. Abrikosov, Sov. Phys. JETP 5 (1957), 1174-1182
  * (2) E. H. Brandt, Rep. Prog. Phys. 58 (1995), 1465-1594
  * (3) E. B. Bogomolny, Sov. J. Nucl. Phys. 24 (1976), 449
  * (4) E. B. Bogomolny and A. I. Vainshtein, Sov. J. Nucl. Phys. 23 (1976), 588-591
  * (5) L. Jacobs and C. Rebbi, Phys. Rev. B 19 (1979), 4486-4494
  * (6) A. J. Beekman and J. Zaanen, Front. Phys. (Beijing) 6 (2011), 357-369 [arXiv:1106.3946 [cond-mat.supr-con]].
  * (7) V. L. Ginzburg and L. D. Landau, Zh. Eksp. Teor. Fiz. 20 (1950), 1064-1082
  * (8) P. G. de Gennes, Superconductivity of Metals and Alloys, (Benjamin, New York, 1966).
  * (9) D. Saint-James, G. Sarma and E. J. Thomas, Type-II Superconductivity, (Benjamin, New York, 1969).
  * (10) G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin and V. M. Vinokur, Rev. Mod. Phys. 66 (1994), 1125-1388
  * (11) M. Tinkham, Introduction to Superconductivity, (Dover Publication, Mineola, New York, 2004).
  * (12) B. Rosenstein and D. Li, Rev. Mod. Phys. 82 (2010), 109-168
  * (13) J. M. Maldacena, Adv. Theor. Math. Phys. 2 (1998), 231-252 [arXiv:hep-th/9711200 [hep-th]].
  * (14) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, Phys. Lett. B 428 (1998), 105-114 [arXiv:hep-th/9802109 [hep-th]].
  * (15) E. Witten, Adv. Theor. Math. Phys. 2 (1998), 253-291 [arXiv:hep-th/9802150 [hep-th]].
  * (16) S. S. Gubser, Phys. Rev. D 78 (2008), 065034 [arXiv:0801.2977 [hep-th]].
  * (17) S. A. Hartnoll, C. P. Herzog and G. T. Horowitz, Phys. Rev. Lett. 101 (2008), 031601 [arXiv:0803.3295 [hep-th]].
  * (18) O. C. Umeh, JHEP 08 (2009), 062 [arXiv:0907.3136 [hep-th]].
  * (19) S. A. Hartnoll, C. P. Herzog and G. T. Horowitz, JHEP 12 (2008), 015 [arXiv:0810.1563 [hep-th]].
  * (20) H. B. Zeng, C. Y. Xia and H. Q. Zhang, JHEP 03 (2021), 136 [arXiv:1912.08332 [hep-th]].
  * (21) E. Nakano and W. Y. Wen, Phys. Rev. D 78 (2008), 046004 [arXiv:0804.3180 [hep-th]].
  * (22) T. Albash and C. V. Johnson, JHEP 09 (2008), 121 [arXiv:0804.3466 [hep-th]].
  * (23) T. Albash and C. V. Johnson, “Phases of Holographic Superconductors in an External Magnetic Field,” [arXiv:0906.0519 [hep-th]].
  * (24) M. Montull, A. Pomarol and P. J. Silva, Phys. Rev. Lett. 103 (2009), 091601 [arXiv:0906.2396 [hep-th]].
  * (25) O. Domenech, M. Montull, A. Pomarol, A. Salvio and P. J. Silva, JHEP 08 (2010), 033 [arXiv:1005.1776 [hep-th]].
  * (26) V. Keranen, E. Keski-Vakkuri, S. Nowling and K. P. Yogendran, Phys. Rev. D 81 (2010), 126012 [arXiv:0912.4280 [hep-th]].
  * (27) T. Albash and C. V. Johnson, Phys. Rev. D 80 (2009), 126009 [arXiv:0906.1795 [hep-th]].
  * (28) K. Maeda, M. Natsuume and T. Okamura, Phys. Rev. D 81 (2010), 026002 [arXiv:0910.4475 [hep-th]].
  * (29) A. Donos, J. P. Gauntlett and C. Pantelidou, JHEP 07 (2020), 095 [arXiv:2001.11510 [hep-th]].
  * (30) A. T. Dorsey, Phys. Rev. B 46 (1992), 8376-8392 [arXiv:cond-mat/9207018 [cond-mat]].
  * (31) Q. Du, M. .D. Gunzburger and J. S. Peterson, Phys. Rev. B 46 (1992) 9027-9034.
  * (32) Ó. J. C. Dias, G. T. Horowitz, N. Iqbal and J. E. Santos, JHEP 04 (2014), 096 [arXiv:1311.3673 [hep-th]].
  * (33) C. Y. Xia, H. B. Zeng, H. Q. Zhang, Z. Y. Nie, Y. Tian and X. Li, Phys. Rev. D 100 (2019) no.6, 061901(R) [arXiv:1904.10925 [hep-th]].
  * (34) X. Li, Y. Tian and H. Zhang, JHEP 02 (2020), 104 [arXiv:1904.05497 [hep-th]].
  * (35) P. Wittmer, C. M. Schmied, T. Gasenzer and C. Ewerz, “Vortex motion quantifies strong dissipation in a holographic superfluid,” [arXiv:2011.12968 [hep-th]].
  * (36) C. Ewerz, A. Samberg and P. Wittmer, “Dynamics of a Vortex Dipole in a Holographic Superfluid,” [arXiv:2012.08716 [hep-th]].
  * (37) A. Srivastav and S. Gangopadhyay, “Novel vortices in a rotating holographic superfluid,” [arXiv:2105.15045 [hep-th]].


# A spatially-resolved large cavity of the J0337 protoplanetary disk in
Perseus

Taichi Uyama Infrared Processing and Analysis Center, California Institute of
Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA  NASA Exoplanet
Science Institute, Pasadena, CA 91125, USA  National Astronomical Observatory
of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Garreth Ruane Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr.,
Pasadena, CA 91109, USA  Kellen Lawson Department of Physics and Astronomy,
University of Oklahoma, Norman, OK 73019, USA  Takayuki Muto  Division of
Liberal Arts, Kogakuin University 2665-1, Nakano-cho, Hachioji-chi, Tokyo,
192-0015, Japan  Charles Beichman NASA Exoplanet Science Institute, Pasadena,
CA 91125, USA  Infrared Processing and Analysis Center, California Institute
of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA  Nienke van
der Marel Leiden Observatory, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands
Physics & Astronomy Department, University of Victoria, 3800 Finnerty Road,
Victoria, BC, V8P 5C2, Canada

###### Abstract

We present Keck/NIRC2 $K_{\rm p}L_{\rm p}$ high-contrast imaging observations
of a J0337 protoplanetary disk. The data discover the spatially-resolved large
cavity, which is the second report among protoplanetary disks in the Perseus
star forming region after the LkH$\alpha$ 330 system. Our data and forward
modeling using RADMC-3D suggests $\sim 80$ au for the cavity radius. There is
discrepancy between J0337’s SED and the modeled SED at $\sim 10\micron$ and
this suggests an unseen inner disk. We also searched for companions around
J0337 but did not detect any companion candidates at separations between
$0\farcs 1$ and $2\farcs 5$. The $L_{\rm p}$-band detection limit corresponds
to $\sim 20M_{\rm Jup}$ at 60 au, $\sim 9-10M_{\rm Jup}$ at 90 au, and $\sim
3M_{\rm Jup}$ at $>120$ au. Compared with other young systems with large
cavities such as PDS 70 and RX J1604, multiple Jovian planets, a single
eccentric Jovian planet, or a massive brown-dwarf at an inner separation could
exist within the cavity.

Protoplanetary Disk, Planet Formation

††journal: AAS Journals

monthyearday []

##  1 Introduction

Perseus is one of the nearby star forming regions ($\sim 300$ pc derived by
Gaia and VLBA observations; Ortiz-León et al., 2018) and has a young age
($\sim 1-3$ Myr for major clusters of IC 348 and NGC 1333; Luhman et al.,
2016). Spectroscopic studies have reported some possible transitional disks in
Perseus - for example, van der Marel et al. (2016) used Spitzer photometry and
IRS data and identified several dozens of transitional disk candidates, some
of which potentially have large cavities ($>$50 au). Transition disk cavities
have been linked to both photoevaporative clearing as part of the evolutionary
disk dissipation process (Alexander et al., 2014) and to the clearing by
massive protoplanets which are already forming in the disk (Lin & Papaloizou,
1979). Previous observational studies suggested Jovian planets in the
transitional disk with cavity $>50$ au (e.g. RX J1604, PDS 70; Dong et al.,
2017; Hashimoto et al., 2012) and eventually VLT/SPHERE high-contrast imaging
observations reported the first convincing protoplanet in the PDS 70 disk
(Keppler et al., 2018). However, the number of convincing protoplanets
embedded in protoplanetary disks is still small and it is important for better
understandings of planet formation mechanisms to detect/characterize more
protoplanets.

2MASS J03370363+3039291 (hereafter J0337) is a member of the Perseus star
forming region, although slightly isolated from IC 348 and NGC 1333, with an
IR excess and van der Marel et al. (2016) suggested a large gap in its disk
from the archival Spitzer catalog. However, other than the gap information a
limited number of stellar/disk parameters are known in this system because
this system has not been prioritized compared with other YSOs in the major
clusters. Here we present Keck/NIRC2 high-contrast imaging observations to
report that our data affirmed the large cavity as predicted in van der Marel
et al. (2016). Section 2 describes our observations and data reduction. In
Section 3 the result of the post-processing is presented. Finally we present
our SED fitting and forward modeling results and discuss the gap opening
mechanisms in the J0337 disk in Section 4.

##  2 Observations and Data Reduction

###  2.1 Observations

J0337 was chosen from the sample of new transition disk candidates with large
cavities from van der Marel et al. (2016), based on SED analysis of hundreds
of Spitzer-selected YSOs. J0337 was selected based on its optical brightness
(PI: Nienke van der Marel). On UT 2016 October 15 we observed J0337 by
Keck/NIRC2 $L_{\rm p}$ band (3.776 µm) combined with the vector vortex mask
and with the total exposure time of 2120 seconds, providing the parallactic
angle change of $91\fdg 33$. After confirming the disk feature (see Section 3)
we conducted follow-up observations for this target using NIRC2 $K_{\rm p}$
band (2.124 µm) on UT 2021 February 2 without a coronagraph mask (PI: Garreth
Ruane, backup target). We achieved 2100 seconds and 38$\fdg$94 for the total
exposure time and the parallactic angle change. Both observations were
conducted under the vertical angle mode for angular differential imaging (ADI;
Marois et al., 2006) to suppress the stellar halo and instrumental speckles.
The typical full-width at half maximum (FWHM) of the point spread function
(PSF) measured 9 pix ($\sim 90$ mas) at $L_{\rm p}$ band and 5 pix ($\sim 50$
mas) at $K_{\rm p}$ band, respectively. We note that no reference-PSF stars
for reference-star differential imaging (RDI; e.g. Ruane et al., 2019) were
observed in both epochs and we apply only ADI for post-processing.

###  2.2 Data Reduction

First we pre-processed (bad pixel correction, flat fielding, sky subtraction,
and image registration) the obtained NIRC2 raw frames (see also Ruane et al.,
2019). We then performed ADI-based post-processing to subtract the stellar
halo by producing the most likely reference PSFs combined with the pyKLIP
algorithms111https://pyklip.readthedocs.io/en/latest/index.html (Wang et al.,
2015) that utilized Karhunen-Loève Image Projection (KLIP; Soummer et al.,
2012). We adopted minrot=10∘ in the pyklip setting, which is not the most
aggressive parameters for the ADI reduction, so that we can avoid a severe
self-subtraction effect that distorts the geometry and attenuates the flux of
the disk features within $\rho\leq 0\farcs 3$ detected in our observations
(see Section 3 for the ADI results).

##  3 Results

In the both epochs we confirmed the same arc-like feature (see Figure 1) with
a signal to noise ratio (SNR) $\sim 4-5$ for its spine, which indicates that
the arc feature is gravitationally bound to J0337. This arc feature is likely
to correspond to forward-scattered light from an inner edge of the outer disk
around J0337. The backward-scatter components are not detected in our
observations. Considering the SNR and the aggressiveness for the ADI reduction
to avoid severe self-subtraction that attenuates the disk feature, we adopted
KL=4 to present our results. However, note that this non-aggressive post-
processing still leaves self-subtraction. We conducted injection test by
burying artificial point sources at a variety of separations and position
angles, which resulted in $>50\%$ self-subtraction particularly in the $K_{\rm
p}$-band data occurs within $0\farcs 3$ because of its small field rotation
angle. The inner dark region indicates a cavity of $\sim 0\farcs 25$ in
radius, which is consistent with the prediction from the SED fitting study
($r_{\rm in}=50^{+30}_{-10}$ au; van der Marel et al., 2016). From the
geometry of the arc feature we roughly estimate an inclination and a position
angle to be $\sim 60^{\circ}$ and $\sim 25^{\circ}$ respectively. The gap size
is estimated by forward modeling (see Section 4.2). Our high-contrast imaging
results show the second large cavity in the Perseus transitional disks after
the LkH$\alpha$ 330 system where a large cavity and a pair of spirals were
reported in its disk (Isella et al., 2013; Akiyama et al., 2016; Uyama et al.,
2018).

---
Figure 1: Keck/NIRC2 observations of J0337 at $K_{\rm p}$ band (left) and
$L_{\rm p}$ band (right). KL=4 is adopted for both images. The central star is
masked by the algorithms. North is up and East is left.

##  4 Discussions

###  4.1 SED Fitting for Stellar Parameters

Due to lack of known spectral type, van der Marel et al. (2016) adopted
$A_{\rm v}=0.0$ but the Bayester19 dust map (Green et al., 2019) suggests
$A_{\rm v}\sim 1.4$ mag around J0337, thus we corrected the photometry data
from published survey data between UV and near-IR wavelengths. Figure 2 shows
J0337’s SED ranging between UV and sub-mm wavelengths (see also Appendix A),
and we note that the previous surveys detected no flux at wavelengths longer
than 200 $\micron$. For SED fitting, we permitted $A_{\rm v}$ values in the
range $1.38-1.46$ (the 1-$\sigma$ range about the nominal value from the
Bayestar19 dust map). Note that we do not account for circumstellar materials
such as envelopes, if any, and in this sense the effective extinction may be
underestimated. We used available GALEX (Martin et al., 2005), TYCHO (Høg et
al., 2000), Gaia (Gaia Collaboration et al., 2016), SDSS (Blanton et al.,
2017), APASS (Henden & Munari, 2014), and 2MASS (Skrutskie et al., 2006)
$JH$-band data points for the stellar characterizations (the best-fit
parameters: $T_{\rm eff}\sim 7800$ K, $\log\ g\sim 5.0$, $L_{\star}\sim 9.3\
L_{\odot}$, and $M_{\star}\sim 1.4\ M_{\odot}$ for $A_{\rm v}=1.38$). SED
retrieval, dereddening, and fitting were carried out using the Virtual
Observatory SED Analyzer (VOSA; Bayo et al., 2008). For SED fitting, we
utilized the BT-Settl-AGSS2009 synthetic stellar spectra (Allard et al., 2012;
Asplund et al., 2009). The VOSA “Chi-square Fit” tool was used to determine
the best-fit parameter values provided above. However, results from other VOSA
SED fitting routines provide comparable solutions (e.g., the “Model Bayes
Analysis” indicates a 1-$\sigma$ confidence interval for $T_{\rm eff}$ in the
range $7444-7800K$). For comparison, we also performed fitting for $A_{\rm
v}=0$ (no extinction; van der Marel et al., 2016); in either case, the derived
best-fit parameters are well located in the 1-5 Myr range of the MIST
evolutionary models (see the left panel of Fig 2, and for details see Dotter,
2016). For further characterizations of the stellar parameters such as
accretion high-dispersion spectroscopy is helpful (Manara et al., 2014; Alcalá
et al., 2017).

We also extrapolated the disk parameters using the updated Gaia-based distance
of 303 pc (Gaia Collaboration et al., 2018). As mentioned in van der Marel et
al. (2016) the estimated disk mass has large uncertainty because the SED
fitting at wavelengths shorter than millimeter can trace only optically-thick
part of the disk. We use an upper limit of the JCMT/SCUBA-2 850 µm
observations, which is converted into $<20M_{\rm Jup}$ with Equation (5) in
van der Marel et al. (2016), for the discussion of the disk evolution
mechanism. In Section 4.2 we also discuss the radiative transfer modeling to
reproduce the outer arc feature. However, we applied ADI technique to remove
stellar halo and to detect the disk scattered light, which attenuates and
distorts the disk feature. Therefore in this study we focus on the detection
of the spatially resolved large cavity. The detailed modeling will be
performed with the follow-up polarimetric differential imaging (PDI; Kuhn et
al., 2001) or RDI data as well as (sub-)millimeter observations of the J0337
disk.

---
Figure 2: (Left) SED of J0337 overlaid with the best-fit stellar photosphere
indicated by the gray hatched region. We used only the data points to the left
of the vertical dashed line to characterize the stellar parameters. The points
to the right of the vertical dashed line are regarded as the IR excess, the
disk parameters for which are estimated in van der Marel et al. (2016). The
inverted triangular point indicates an upper limit. The archival photometry is
summarized in Appendix A. (Right) Comparison of MIST evolutionary tracks
(solid lines) and isochrones (dashed lines) with the best-fit parameters of
J0337’s SED, assuming either the nominal $A_{\rm V}$ value or no extinction.
These parameters are consistent with a 1–5 Myr star, in agreement with the age
estimate of the Perseus star forming region.

###  4.2 Forward Modeling

To validate the cavity size we forward-modeled the disk feature to the NIRC2
data using RADMC-3D radiative transfer code (Dullemond et al., 2012) and
pyklip forward modeling modules. As mentioned in Section 4.1, however, it is
difficult to estimate the detailed disk parameters with the current SED data
points and we aim at reproducing the NIRC2 image and the MIR excess of the SED
in the forward modeling processes. For simplicity we fix the
inclination$=60^{\circ}$ and position angle$=-25^{\circ}$ in the modeling. We
found that the inner wall of the outer disk needs to have low surface density
so that the forward scattering is brighter than the backward scattering with
inclination$\sim 60^{\circ}$, which can reproduce the NIRC2 results, on the
other hand the outer disk needs to have a high enough dust mass (high surface
density) to reproduce the MIR excess. A power-law density profile
($\Sigma_{\rm dust}\propto r^{-0.5}$, where $\Sigma_{\rm dust}$ and $r$
correspond to the dust surface density and separation from the star
respectively) can reproduce either of these two characteristics but cannot
reproduce both simultaneously. Therefore, for reproducing both characteristics
we utilized a Gaussian profile

$\displaystyle\Sigma_{\rm dust}(r)=\Sigma_{\rm center}*\exp{(-(r-r_{\rm
center})^{2}/2\sigma_{r}^{2})},$ (1)

where $\Sigma_{\rm center}=5\times 10^{-3}$ g cm-2, $r_{\rm center}$= 200 au,
and $\sigma_{r}$= 50 au; note that these parameters are degenerate). In this
paper, we fix $\Sigma_{\rm center}$, $r_{\rm center}$, and $\sigma_{r}$ since
near-IR does not in general trace dust surface density. The dust surface
density distribution will be better constrained by spatially resolved sub-mm
observations with e.g., ALMA.

We also adopted a dust settling factor $f_{\rm set}=2$ to match the
observational results and a scale height is calculated as $\frac{\sqrt{kT_{\rm
dust}/(2.34m_{\rm H})}}{f_{\rm set}\sqrt{GM_{\star}/r^{3}}}$, where $k,T_{\rm
dust},m_{\rm H},G$ are the Boltzmann constant, dust temperature, mass of
hydrogen, and gravitational constant, respectively. We used thermal Monte
Carlo computations of RADMC-3D (mctherm) to calculate the dust temperature
derived from the adopted dust density profile and the stellar parameters. We
then utilized the Mie scattering codes (Bohren & Huffman, 1983) with the
Gaussian dust-size distribution (centered around 10 $\micron$, 5% width) and
the optical constant of amorphous silicate222http://www.astro.uni-
jena.de/Laboratory/OCDB/data/silicate/amorph/pyrmg70.lnk to calculate dust
opacity. These assumptions provide a good match to the catalog values around
MIR and we tested the cavity size by defining a cutoff of the modeled disk
($\Sigma_{\rm dust}(r)=0$ at $r<r_{\rm cutoff}$ in Equation 1) as the cavity
radius and changing it from 40 au to 100 au by a 20-au step. We injected the
modeled disk to the NIRC2 data at a position angle of -115∘ (rotated by 90∘
clockwise) and then rereduced them. The 80-au cavity model seems to best match
the NIRC2 results while the SEDs are not largely affected within the parameter
range of the cavity radius we explored (see Figures 3 and 4 and Appendix B).
Note that because of the imperfect AO correction and the presence of the
actual disk feature, the injected disk feature is affected to some extent in
the post-processed image. The $L_{\rm p}$-band forward-modeled images are less
clear than the $K_{\rm p}$-band images as the $L_{\rm p}$-band observations
have the larger PSF (see also Section 2.1). Particularly the southern part of
the injected disk is attenuated by self-subtraction of the actual disk
features. In Table 1 we summarize the stellar and disk parameters we adopted
in this study. Since the system is yet relatively unexplored and we only have
images in near-IR wavelengths, it is hard determine all the model parameters.
We therefore address the uncertainties of the model parameters briefly by
changing stellar luminosity, whose uncertainty is expected to be a factor of a
few due to the uncertainty of extinction towards the star. The change of the
stellar luminosity affects the disk temperature. As shown in Appendix C, we
see that the modeled SED is affected by the stellar luminosity while the
geometry of the modeled near-IR image is not. Therefore, we consider that the
cavity radius of $\sim 80$ au is relatively robust despite the uncertainty of
stellar luminosity.

---
Figure 3: The adopted model assuming the cavity size of 80 au. (Left) The
modeled SED in this study compared with the catalog values - we include only
the stellar component with $T_{\rm eff}=7800$ K blackbody and the outer disk
component (see text). The discrepancy at $\sim 10\micron$ suggests an unseen
inner disk. (Middle) The modeled disk image at $K_{\rm p}$ band overlaid with
the NIRC2 contours. (Right) Same as the middle panel for $L_{\rm p}$ band.
Note that the modeled disk images are smoothed by Gaussian with the measured
FWHM of each NIRC2 observation.
---
Figure 4: Post-processed images at $K_{\rm p}$ (left) and $L_{\rm p}$ band
(right) after injecting the 80au-cavity model (Figure 3) at position angle of
-115∘.

We also note that our modeled SED cannot reproduce the $\sim 10\micron$ excess
and the discrepancy indicates presence of an unseen inner disk. As the NIRC2
data did not confirm any inner disk features the inner disk may be located
within or close to the inner working angle of the NIRC2 observations ($\sim
0\farcs 1=30\ {\rm au}$), which is indeed well outside the typical radii of
inner disks (Francis & van der Marel, 2020). Future interferometric
observations or follow-up PDI/RDI observations will be able to discuss the
inner disk. Particularly the follow-up PDI/RDI may detect a shadowing effect
by the inner disk (e.g. Bohn et al., 2021) and this should be taken into
account for a better model of the J0337 disk.

Table 1: Summary of the adopted parameters parameter | value
---|---
Stellar parameters
$T_{\rm eff}$ | 7800 K
$\log g$ | 5.0
$L_{\star}$ | 9.3 $L_{\odot}$
$M_{\star}$ | 1.4 $M_{\odot}$
age | 1-5 Myr
$A_{\rm V}$ | 1.38
Disk parameters
inclination | 60∘
position angle | -25∘
dust distribution | Gaussian profile (Eqn. 1)
cavity radius | 80 au

###  4.3 Gap Opening Mechanism

Alexander et al. (2006) simulated the disk evolution with photoevaporation and
predicted that the disk with a photoevaporation-induced cavity has the order
of $M_{\rm disk}\sim 10^{-4}\ M_{\star}$. With the updated stellar parameters
and the upper limit of the disk mass (see Section 4.1), the disk mass is
$<1.4\times 10^{-2}\ M_{\star}$. If the actual disk mass is smaller than the
upper limit by an order of magnitude photoevaporation can be responsible for
the cavity, but it is unlikely that photoevaporation is the only mechanism
that opens such a large cavity (e.g. Espaillat et al., 2010). Accretion rate
is another index to discuss the mass loss ratio by photoevaporation (e.g. Owen
et al., 2012; Ercolano & Pascucci, 2017), but J0337’s accretion has not been
studied. Deep radio-wavelengths explorations and high-dispersion spectroscopy
for this system will help to determine the disk mass and to further
investigate the photoevaporation scenario.

Planet formation is also one of the most plausible scenarios given that the
cavity size is similar to the PDS 70 system (Keppler et al., 2018). However,
our observations did not detect any convincing companion candidates at
separations between $0\farcs 1$ and $2\farcs 5$. Figure 6 shows the $L_{\rm
p}$-band ADI result with KL=10 to search for companions. Figure 6 shows
5$\sigma$ detection limit of Figure 6 with the expected contrast of a
substellar-mass object assuming 1 Myr and COND03 model (Baraffe et al., 2003).
We divided the FoV into a number of annular regions and defined the noise as
the standard deviation within each annular region, where we masked the disk
feature and used 5$\sigma$ clipping to mitigate the effect of the disk feature
on the detection limit. We also corrected the self-subtraction effect by
injecting fake sources at different position angle and separations and
calculating the flux-attenuation ratio. Out detection limits could constrain
the presence of a massive brown-dwarf companion in the cavity region ($\sim
20\ M_{\rm Jup}$ at 60 au and $\sim 9-10\ M_{\rm Jup}$ at 90 au) and that of a
$\sim 3\ M_{\rm Jup}$ protoplanet outside the cavity ($>120$ au). Compared
with other high-contrast imaging and simulation studies that targeted large-
cavity transitional disks such as PDS 70 or RX J1604 (Haffert et al., 2019;
Muley et al., 2019; Canovas et al., 2017) there could be multiple Jovian
planets less than $10\ M_{\rm Jup}$, a single eccentric Jovian planet, or a
brown-dwarf companion close to the central star that we cannot resolve in our
observations.

Figure 5: Same NIRC2 $L_{\rm p}$-band result (see Figure 1) with the larger
FoV and KL=10 to search for outer companion candidates. Figure 6: 5$\sigma$
detection limit of the $L_{\rm p}$-band observations made from the KL=10
result (see Figure 6).
---

##  5 Summary

YSOs in Perseus are typically more distant than other nearby star forming
regions such as Taurus, Lupus, and Ophiuchus. Furthermore Perseus is located
at northern in the celestial field and Perseus has been less prioritized by
ALMA or adaptive optics observations. In this study we presented Keck/NIRC2
$K_{\rm p}L_{\rm p}$ high-contrast imaging observations of the J0337
protoplanetary disk, one of the transitional disk candidates with a large
cavity suggested in van der Marel et al. (2016). After ADI reduction with
pyklip we detected forward scattering from an inner wall of the outer disk and
affirmed the large cavity. As the ADI reduction distorts the geometry of the
disk, we investigated the cavity size with forward-modeled disks that are made
from RADMC-3D radiative transfer codes as well as aiming to reproduce its SED.
Our data and forward modeling suggest $\sim 80$ au for the cavity radius,
which is consistent with the prediction of van der Marel et al. (2016).
However there is discrepancy at $\sim 10\micron$ in J0337’s SED and the
modeled SED suggesting the presence of an unseen inner disk. We also searched
for companions around J0337 but did not detect any convincing companions
within $2\farcs 5$ and the $L_{\rm p}$-band detection limit could set a
constraint on the mass of potential companions to $\sim 20M_{\rm Jup}$ at 60
au, $\sim 9-10M_{\rm Jup}$ at 90 au, and $\sim 3M_{\rm Jup}$ at $>120$ au
assuming 1 Myr and COND03 model. Compared with other protoplanetary disk
systems with large cavities such as PDS 70 and RX J1604, a massive brown-dwarf
companion at an inner separation or multiple Jovian planets fainter than our
detection limit within the cavity may exist inside the cavity.

The detailed modeling, including the inner disk component, will help better
understand the J0337 system while follow-up observations such as spectroscopic
observation of J0337 for the stellar characterization, PDI/RDI to investigate
the scattered light without self-subtraction, obtaining radio flux, and ALMA
high-angular resolution observation are essential. Finally our discovery is
the second spatially resolved cavity in the Perseus protoplanetary disks and
there remain more transitional disks suggested in SED-based studies such as
van der Marel et al. (2016), which promotes further high angular resolution
observations at this region.

The authors would like to thank the anonymous referees for their constructive
comments and suggestions to improve the quality of the paper. We thank Akimasa
Kataoka for helpful comments to develop the RADMC-3D modeling. T.U. is
supported by Grant-in-Aid for Japan Society for the Promotion of Science
(JSPS) Fellows and JSPS KAKENHI Grant No. JP21J01220. T. M. is supported by
JSPS KAKANHI Grant Nos. 17H01103, 18H05441 and 19K03932. The data presented
herein were obtained at the W. M. Keck Observatory, which is operated as a
scientific partnership among the California Institute of Technology, the
University of California and the National Aeronautics and Space
Administration. The Observatory was made possible by the generous financial
support of the W. M. Keck Foundation. The authors wish to recognize and
acknowledge the very significant cultural role and reverence that the summit
of Maunakea has always had within the indigenous Hawaiian community. We are
most fortunate to have the opportunity to conduct observations from this
mountain. This research has made use of the Keck Observatory Archive (KOA),
which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science
Institute (NExScI), under contract with the National Aeronautics and Space
Administration. This publication makes use of VOSA, developed under the
Spanish Virtual Observatory project supported by the Spanish MINECO through
grant AyA2017-84089. VOSA has been partially updated by using funding from the
European Union’s Horizon 2020 Research and Innovation Programme, under Grant
Agreement nº 776403 (EXOPLANETS-A)

##  Appendix A J0337 SED

Table 2 summarizes archival photometry of J0337 (see also Figure 2). We used
these catalog values no longer than 1.7$\micron$ for the SED fitting to
estimate the stellar parameters and referred to the MIR photometry for the
forward modeling.

Table 2: Archival photometry of J0337 wavelength [µm] | flux [erg s-1 cm-2 µm-1] | flux error [erg s-1 cm-2 µm-1] | remarks
---|---|---|---
0.155 | 2.2e-10 | 3.97e-11 | GALEX (Martin et al., 2005)
0.23 | 4.23e-09 | 6.46e-11 | GALEX (Martin et al., 2005)
0.428 | 5.75e-09 | 6.63e-10 | TYCHO (Høg et al., 2000)
0.43 | 6.60e-09 | 2.61e-10 | APASS (Henden & Munari, 2014)
0.467 | 6.00e-09 | 2.65e-10 | SDSS (Blanton et al., 2017)
0.534 | 4.43e-09 | 4.24e-10 | TYCHO (Høg et al., 2000)
0.539 | 4.82e-09 | 4.22e-10 | APASS (Henden & Munari, 2014)
0.582 | 3.11e-09 | 8.17e-12 | Gaia (Gaia Collaboration et al., 2016)
0.614 | 3.69e-09 | 1.87e-10 | SDSS (Blanton et al., 2017)
0.746 | 2.36e-09 | 1.67e-10 | SDSS (Blanton et al., 2017)
1.235 | 5.83e-10 | 1.23e-11 | 2MASS (Skrutskie et al., 2006)
1.662 | 2.44e-10 | 4.71e-12 | 2MASS (Skrutskie et al., 2006)
2.159 | 1.11e-10 | 1.74e-12 | 2MASS (Skrutskie et al., 2006)
3.353 | 4.44e-11 | 9.41e-13 | WISE (Wright et al., 2010)
4.603 | 1.91e-11 | 3.51e-13 | WISE (Wright et al., 2010)
8.228 | 1.22e-11 | 5.39e-13 | AKARI (Ishihara et al., 2010)
10.146 | 8.98e-12 | 1.08e-12 | IRAS (Neugebauer et al., 1984)
11.561 | 4.21e-12 | 6.2e-14 | WISE (Wright et al., 2010)
17.609 | 4.89e-12 | 2.26e-13 | AKARI (Ishihara et al., 2010)
21.727 | 7.16e-12 | 5.85e-13 | IRAS (Neugebauer et al., 1984)
22.088 | 4.83e-12 | 1.11e-13 | WISE (Wright et al., 2010)
23.21 | 5.42e-12 | 5.04e-13 | Spitzer (Evans et al., 2003)
51.989 | 6.22e-12 | 5.24e-13 | IRAS (Neugebauer et al., 1984)
62.951 | 3.26e-12 | 3.6e-14 | AKARI (Ishihara et al., 2010)
68.445 | 2.57e-12 | 2.48e-13 | Spitzer (Evans et al., 2003)
76.904 | 2.46e-12 | 1.62e-13 | AKARI (Ishihara et al., 2010)
95.297 | 1.62e-12 | 7.56e-13 | IRAS (Neugebauer et al., 1984)
140.856 | 4.62e-13 | 5.26e-14 | AKARI (Ishihara et al., 2010)
857.914 | 2.27e-16 | 1.17e-16 | JCMT upper limit (van der Marel et al., 2016)

##  Appendix B Forward Modeling with Different Cavity Sizes

Figures 7, 8, and 9 compare the model images and the forward-modeled results
assuming 40, 60, and 100 au for the cavity size (defined as a cutoff in the
model) in radius respectively. The Gaussian radial profile for the outer disk
is fixed in the modeling. Regarding the scattering profile we multiplied the
$K_{\rm p}$ and $L_{\rm p}$ modeled disks by 0.8 and 2 to correct the
difference between the simulated blackbody with $T_{\rm eff}=7800$ K and the
catalog values. Considering the geometry of the disk scattered light, the
80au-cavity model best matches the NIRC2 results. In the $K_{\rm p}$-band
forward-modeled images, the 40au- and 60au-cavity models are reproduced as an
almost straight line because the ADI reduction distorted the modeled disk,
while the NIRC2 data show the arc feature as presented in Section 3. The
100au-cavity model shows the larger ring feature than the NIRC2 data and
suggests the bright backward-scattering feature that is not seen in fact.

---
Figure 7: The modeled disk images assuming the cavity size of 40 au. From left
to right: modeled disk at $K_{\rm p}$ band overlaid with the NIRC2 result
contours (black), post-processed NIRC2 image at $K_{\rm p}$ band after
injecting the modeled disk (forward-modeled result), modeled disk at $L_{\rm
p}$ band, and forward-modeled result at $L_{\rm p}$ band.
---
Figure 8: Same as Figure 7 for the case of 60au-cavity.
---
Figure 9: Same as Figures 7, 8 for the case of 100au-cavity.

##  Appendix C RADMC-3D Modeling with Different Dust Temperature

The uncertainty on luminosity from the VOSA fitting is 0.18 but as mentioned
in Section 4.1 the uncertainty on the extinction estimation could affect
characterizing the stellar parameters. To test the effect of the dust
temperature on the RADMC-3D modeling results, we varied a larger range of the
stellar luminosity from 3.1 $L_{\odot}$ ($L_{\star}/3$) to 27.9 $L_{\odot}$
($3L_{\star}$) by changing the stellar temperature in the RADMC-3D settings
from 5927 ($7800/3^{0.25}$) K to 10265 ($7800\times 3^{0.25}$) K while we
fixed the cavity size to 80 au. The dust temperature was accordingly
calculated by the ’mctherm’ command in RADMC-3D. Figure 10 shows the modeled
SEDs and Figure 11 compares the modeled images with the above assumptions. The
SED and surface brightness of the disk are variable with the stellar
luminosity while the geometry of the disk feature does not vary within this
luminosity range.

Figure 10: Comparison of the modeled SEDs assuming different temperature of
the central star and the 80-au cavity. The SED with $T_{\star}=7800$ K is the
same as Figure 3.
---
Figure 11: Comparison of the modeled images assuming different temperature of
the central star and the 80-au cavity. From left to right: modeled disk at
$K_{\rm p}$ band with $T_{\star}=5926$ K, modeled disk at $L_{\rm p}$ band
with $T_{\star}=5926$, and modeled disk at $K_{\rm p}$ and $L_{\rm p}$ bands
with $T_{\star}=10265$ K.

## References

  * Akiyama et al. (2016) Akiyama, E., Hashimoto, J., Liu, H. B., et al. 2016, AJ, 152, 222, doi: 10.3847/1538-3881/152/6/222
  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929
  * Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, Protostars and Planets VI, 475, doi: 10.2458/azu_uapress_9780816531240-ch021
  * Alexander et al. (2006) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229, doi: 10.1111/j.1365-2966.2006.10294.x
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765, doi: 10.1098/rsta.2011.0269
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Baraffe et al. (2003) Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701, doi: 10.1051/0004-6361:20030252
  * Bayo et al. (2008) Bayo, A., Rodrigo, C., Barrado Y Navascués, D., et al. 2008, A&A, 492, 277, doi: 10.1051/0004-6361:200810395
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/aa7567
  * Bohn et al. (2021) Bohn, A. J., Benisty, M., Perraut, K., et al. 2021, arXiv e-prints, arXiv:2112.00123.  https://arxiv.org/abs/2112.00123
  * Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles
  * Canovas et al. (2017) Canovas, H., Hardy, A., Zurlo, A., et al. 2017, A&A, 598, A43, doi: 10.1051/0004-6361/201629145
  * Dong et al. (2017) Dong, R., van der Marel, N., Hashimoto, J., et al. 2017, ApJ, 836, 201, doi: 10.3847/1538-4357/aa5abf
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool.  http://ascl.net/1202.015
  * Ercolano & Pascucci (2017) Ercolano, B., & Pascucci, I. 2017, Royal Society Open Science, 4, 170114, doi: 10.1098/rsos.170114
  * Espaillat et al. (2010) Espaillat, C., D’Alessio, P., Hernández, J., et al. 2010, ApJ, 717, 441, doi: 10.1088/0004-637X/717/1/441
  * Evans et al. (2003) Evans, Neal J., I., Allen, L. E., Blake, G. A., et al. 2003, PASP, 115, 965, doi: 10.1086/376697
  * Francis & van der Marel (2020) Francis, L., & van der Marel, N. 2020, ApJ, 892, 111, doi: 10.3847/1538-4357/ab7b63
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749, doi: 10.1038/s41550-019-0780-5
  * Hashimoto et al. (2012) Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJ, 758, L19, doi: 10.1088/2041-8205/758/1/L19
  * Henden & Munari (2014) Henden, A., & Munari, U. 2014, Contributions of the Astronomical Observatory Skalnate Pleso, 43, 518
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Isella et al. (2013) Isella, A., Pérez, L. M., Carpenter, J. M., et al. 2013, ApJ, 775, 30, doi: 10.1088/0004-637X/775/1/30
  * Ishihara et al. (2010) Ishihara, D., Onaka, T., Kataza, H., et al. 2010, A&A, 514, A1, doi: 10.1051/0004-6361/200913811
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44, doi: 10.1051/0004-6361/201832957
  * Kuhn et al. (2001) Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJ, 553, L189, doi: 10.1086/320686
  * Lin & Papaloizou (1979) Lin, D. N. C., & Papaloizou, J. 1979, MNRAS, 188, 191
  * Luhman et al. (2016) Luhman, K. L., Esplin, T. L., & Loutrel, N. P. 2016, ApJ, 827, 52
  * Manara et al. (2014) Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18, doi: 10.1051/0004-6361/201323318
  * Marois et al. (2006) Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556, doi: 10.1086/500401
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1, doi: 10.1086/426387
  * Muley et al. (2019) Muley, D., Fung, J., & van der Marel, N. 2019, ApJ, 879, L2, doi: 10.3847/2041-8213/ab24d0
  * Neugebauer et al. (1984) Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJ, 278, L1, doi: 10.1086/184209
  * Ortiz-León et al. (2018) Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880, doi: 10.1111/j.1365-2966.2011.20337.x
  * Ruane et al. (2019) Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118, doi: 10.3847/1538-3881/aafee2
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Soummer et al. (2012) Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28, doi: 10.1088/2041-8205/755/2/L28
  * Uyama et al. (2018) Uyama, T., Hashimoto, J., Muto, T., et al. 2018, AJ, 156, 63, doi: 10.3847/1538-3881/aacbd1
  * van der Marel et al. (2016) van der Marel, N., Verhaar, B. W., van Terwisga, S., et al. 2016, A&A, 592, A126, doi: 10.1051/0004-6361/201628075
  * Wang et al. (2015) Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, in Astrophysics Source Code Library, ascl:1506.001.  http://ascl.net/1506.001
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868