# Causal Analysis on the Anchor Store Effect in a Location-based Social
Network

Anish K. Vallapuram1, Young D. Kwon2, Lik-Hang Lee3, Fengli Xu4, and Pan Hui1
Email<EMAIL_ADDRESS>1Hong Kong University of Science and
Technology, Hong Kong SAR  2University of Cambridge, United Kingdom  3Korea
Advanced Institute of Science and Technology, South Korea  4University of
Chicago, United States

###### Abstract

A particular phenomenon of interest in Retail Economics is the spillover
effect of anchor stores (specific stores with a reputable brand) to non-anchor
stores in terms of customer traffic. Prior works in this area rely on small
and survey-based datasets that are often confidential or expensive to collect
on a large scale. Also, very few works study the underlying causal mechanisms
between factors that underpin the spillover effect. In this work, we analyze
the causal relationship between anchor stores and customer traffic to non-
anchor stores and employ a propensity score matching framework to investigate
this effect more efficiently. First of all, to demonstrate the effect, we
leverage open and mobile data from London Datastore and Location-Based Social
Networks (LBSNs) such as Foursquare. We then perform a large-scale empirical
analysis of customer visit patterns from anchor stores to non-anchor stores
(e.g., non-chain restaurants) located in the Greater London area as a case
study. By studying over 600 neighbourhoods in the Greater London area, we find
that anchor stores cause a 14.2-26.5% increase in customer traffic for the
non-anchor stores reinforcing the established economic theory Moreover, we
evaluate the efficiency of our methodology by studying the confounder balance,
dose difference and performance of the matching framework on synthetic data.
Through this work, we point decision-makers in the retail industry to a more
systematic approach to estimate the anchor store effect and pave the way for
further research to discover more complex causal relationships underlying this
effect with open data.

††publicationid: pubid:  IEEE/ACM ASONAM 2022, November 10-13, 2022
http://dx.doi.org/10.1145/XXXXXXX.XXXXXXX 978-1-6654-5661-6/22/$31.00 © 2022
IEEE

##  I Introduction

Retail store co-location is a widely studied problem in the field of Economics
[1, 2, 3, 4]. A particular phenomenon of interest is the _anchor store
effect_. The idea of an anchor store originates from large department stores
which cover a wide range of products which causes a spillover of shopper
footfall within its vicinity (i.e., to non-anchor stores). This can occur due
to multitudinous factors including but not limited to the brand name, consumer
awareness of the product prices and their reduced commute costs per product
category. Anchor stores are advantageous to owners of small businesses and
malls alike as they can leverage the footfall from the anchor stores’ brand
name. Motivated by this, mall owners offer discounted rents to anchor stores
for being their tenants [5].

Previous studies mainly focus on the effect of large anchor stores in malls on
the sales and customer traffic of small non-anchor tenants [3]. Several works
[6, 7] have investigated correlations between the rents paid by both anchor
and non-anchor tenants to the sales. Some other works [8, 9] have conducted
surveys on the prospect that customers in a neighbourhood will visit a mall
given the brand reputation of its anchor tenants. However, these works only
consider mall-related variables which are insufficient for the study as later
empirical studies [10, 11, 12] demonstrate the effect of the entrance of
anchor stores on the retail trade at the town level. Furthermore, Daunfeldt et
al. [13] expostulate that these works employ descriptive statistics to
quantify the effects and do not draw causal conclusions. Still, these works
examine the impact of anchor stores on retail sales, profit margins and
employment. The data for these variables is either under non-disclosure
agreements or cost-ineffective to acquire at a large scale. This limits the
studies to a specific retail brand, ignoring the other anchor stores in the
vicinity.

In this paper, we instead leverage publicly available user mobility data
generated by the location-based social networks (LBSNs), Foursquare, and
census data London Datastore which covers observational data for 625
neighbourhoods across Greater London over three years. We then estimate the
causal impact (i.e., change in customer traffic) of the presence of anchor
stores (e.g. Tesco, Sainsbury’s, Waitrose, Lidl) on non-anchor stores located
in the Greater London area as a case study (Figure 1). In our initial causal
analysis, we focus on non-chain restaurants as a representative example of
non-anchor stores to narrow down the scope of this study. Yet, our method is
generalisable to any other type of non-anchor store. Furthermore, we go beyond
the mall-related variables and consider several neighbourhood characteristics
that significantly impact customer traffic to retail stores [7]. Additionally,
prior works either ignore these variables or fail to consider their
interdependence leading to biased conclusions on the effect. In contrast, our
work mitigates the aforementioned issue by employing a causality analysis
framework that handles the interdependence between variables before performing
the inference and drawing conclusions. To the best of our knowledge, this is
the first work to leverage data from LBSNs and employ causal inference to
demonstrate the effect of anchor stores on increased customer traffic to
neighbouring non-anchor stores.

Overall, the major contributions of this paper are:

  * •

To the best of our knowledge, this is the first work to employ observational
data from an LBSN to perform a large-scale causal study on this effect with
over 1,583 anchor stores and 7,694 non-anchor restaurants in the Greater
London area by observing customer traffic over three years.

  * •

Based on a propensity score matching framework, we study the effectiveness of
the anchor store effect by considering neighbourhood characteristics such as
population density, disposable income, and transport accessibility to further
mitigate their confounding effect on customer traffic to anchor stores and
retail restaurants.

  * •

Our results consistently demonstrate a 14.2-26.5% increase in customer traffic
in non-anchor stores with more neighbouring anchor stores, reinforcing the
well-known economic concept. We further evaluate the significance of these
results through several metrics and performance on synthetic data with known
causal effects.

  * •

Based on the findings of this work on the causal relation between anchor
stores and retail restaurants, we elaborate on how our findings are useful to
various beneficiaries.

Figure 1: Our problem formulation – studying the anchor store effects on non-
chain restaurants. The number of anchor stores per restaurant (left) and the
number of tips per restaurant (right) are shown for each neighbourhood in the
Greater London area. Data in 2012 are depicted, with the darker colour
representing a higher count of stores/tips per restaurant or vice versa.

##  II Related Work

We will first review previous works on anchor store effects and then delve
into the background of causality analysis.

Anchor Store Effect. Several works have investigated the effect of anchor
stores in malls on the sales and customer traffic to non-anchor tenants [3].
[6, 7] study this impact by finding the correlation between the rents paid by
anchor stores and non-anchor stores and [4] perform regression analysis on the
mall visitor traffic with respect to several independent variables such as the
total area of the commercial centre, age, number of stores, number of anchor
stores and their total area. [7] suggests that variables beyond the walls of
the mall also affect customer traffic. Consequently, several empirical studies
[10, 11, 12] demonstrated the effect of the entrance of anchor stores on
retail trade at the town level. For example, Davidson et al. [10] concluded
that total revenues increase in towns by 41% in Wal-Mart entry towns versus
28% in towns which Wal-Mart did not enter. Yet, these works do not draw causal
conclusions as they quantify the effect using descriptive statistics. Hence,
Daunfeldt et al. [13] employ a propensity score matching framework [14] to
study the effect of the entrance of IKEA stores in three municipalities in
Sweden and concluded a 7% increase in revenues to incumbent retailers.

However, these works study the impact of anchor stores on variables such as
revenues, profit margins, wages, etc. which poses two key limitations: 1) the
data is not available for analytics due to non-disclosure agreements, and 2)
The data collection is laborious and cost-ineffective at a large scale. These
works limit themselves to a specific retail brand, ignoring the other anchor
stores in the vicinity. Thus, in this work, we study the impact of anchor
stores on customer traffic by employing observational data from an LBSN which
allows us to study the impact of 1,583 anchor stores from several retail
brands on 7,694 restaurants in the Greater London area.

Causality Analysis. To ensure our findings’ robustness and generalisability,
we study the anchor store effect via causality analysis, a long-standing
research area that lies at the heart of scientific discovery [15]. Randomised
controlled trial (RCT) is a widely accepted approach for performing this
analysis, especially in medicine and biology [16]. However, it requires
randomly assigning experiment subjects to treated and controlled groups in
customised experiments, which is often expensive and infeasible. Several
alternative frameworks have addressed this by estimating the causal
relationships from observational data, which fall into three categories: 1)
Given the underlying causal diagram of observational data, causal relation can
be estimated with statistical methods, e.g., instrumental variables analysis
[17, 18]. However, such a causal diagram requires extensive expert knowledge
and is difficult to validate; 2) Under the assumption that we observe all the
confounding variables, simple regression analysis can examine the causality in
observational data [19], but these approaches are sensitive to the selected
functional forms of the confounders. Even if we could observe all confounders,
we still face the common support problem where the covariate space must align
across all the subjects. 3) Another angle is to estimate the result of RCT by
matching each treated subject with a controlled subject that has similar
confounding variables [20]. In this case, the complexity of the common support
problem is limited to the subjects in the matched pair.

In this paper, thus design a non-randomised study that estimates the causal
effect of an anchor store with a graphical propensity score matching
algorithm, which falls in the third category mentioned above.

##  III Methodology

Formally, the goal of this work is to study the effect of the presence of
anchor stores within the neighbourhood of non-anchor stores on their customer
visit patterns. In a causal study, the units are the basic experimental
objects belonging to a larger collection of the population which in our case
are neighbourhoods in a large city. A highly relevant instance is the Greater
London area where the neighbourhoods are analogous to its 625
wards111neighbourhoods, units and wards will now be used interchangeably. The
presence of the anchor stores is the treatment denoted by $Z$, and customer
visit patterns to non-anchor stores, such as non-chain restaurants, is the
outcome denoted by $Y$.

We perform causality analysis by employing Rubin’s framework of potential
outcomes [15] which operates as follows. Given a set of units $U$ and a set of
binary treatment levels $Z\in\\{0,1\\}$, we are interested in comparing the
outcome of applying each treatment level $Y_{1}(u)$ and $Y_{0}(u)$ for each
unit $u$. The average treatment effect (ATE) is then $E(Y_{1}(u)-Y_{0}(u))$.
Yet, the rudimentary problem of causal analysis is that both the outcomes for
a unit cannot be observed, especially in our context. Therefore, a widely
sought-after alternative [21, 22] is a matching-based approach that determines
the effect of the treatment $Z$ on the outcome $Y$ by comparing pairs of units
with highly similar confounder variables or simply confounders $X$. These
variables describe certain fundamental characteristics of the units that have
a causal effect on the treatment and outcome variables themselves. For
example, customer visit patterns to non-chain restaurants (outcome) can be
affected by not only the presence of the anchor stores (treatment) but also by
the population density of the neighbourhood (confounder). Ignoring the effect
of confounders is a fundamental fallacy in causality analysis. Once the
confounders are known, a similarity of distance metric for two units can be
defined based on their confounders. The final set of matching $M$ contains
highly similar pairs of units $(u_{t},u_{c})$ where the treatment values of
the pair are $Z(u_{t})=1$ and $Z(u_{c})=0$ and no unit appears in more than
one pair. The average treatment effect (ATE) can be computed as

$ATE=\frac{1}{|M|}\sum_{(u_{t},u_{c})\in M}Y(u_{t})-Y(u_{c})$ (1)

The remainder of this section addresses how we adapt the matching framework to
our problem based on one remark. The causal framework has been described with
binary treatment levels, i.e. $Z\in(0,1)$. However, the treatment variable
takes a range of values in observational data for our problem and thus binned
into treatment levels or doses as shown in Table I.

###  III-A Matching

This section details a graphical approach to match units based on confounder
similarity. We first construct a fully connected graph $G=(V,E)$ where the
vertices are units. Each edge $(u_{i},u_{j},w_{ij})\in E$ has a weight
($w_{ij}$) determined by the distance metric based on the confounders of the
units ($u_{i},u_{j}$). An optimal matching $M\subset E$ is achieved by
optimising:

$\operatorname*{argmin}_{M}\frac{1}{|M|}\sum_{(u_{i},u_{j},w_{ij})\in
M}w_{ij}$ (2)

where the maximum matching must be determined by the lowest sum of edge
weights. A matching $M$ is called a maximum matching if it includes the most
edges from the graph, i.e. for every matching $M^{\prime}$ of $G$,
$|M|\geq|M^{\prime}|$. In a causal experiment, outcomes are only compared
between units from different treatment groups [20]. This can be achieved by
inducing a new graph $G^{\prime}=(V,E^{\prime})$ from $G$ by simply dropping
all edges in $E$ where the connected units have the same treatment level, i.e.
$E^{\prime}=\\{(u_{i},u_{j},w_{ij})\in E:z_{i}>z_{j}$.

###  III-B Propensity Score

An appropriate choice of distance metric is paramount for the quality of the
matching. While euclidean distances appear as obvious choices, they do not
scale well with the dimensionality of the confounders and there will be non-
uniform confounder distributions within the disjoint vertex sets which not
ideal. An alternative metric called propensity score resolves this [14], which
is defined as the probability that a unit receives a particular treatment
level given the confounders. The propensity score function outputs the
probability distribution for each treatment level, hence the dimensionality is
more controlled compared to the dimensionality of the confounders. Joffe and
Rosenbaum [23] show that the propensity score must satisfy the balancing
property
$P(X=X(u_{i})|e(X(u_{i})),Z=z)=P(X=X(u_{i})|e(X(u_{i})),Z=z^{\prime})$ for all
combinations of treatment levels $(z,z^{\prime})$ where $e$ is the propensity
score function. This is certainly the case when the propensity score is
described using the McCullagh’s ordered logit model [24]

$log\left(\frac{P(z_{i}>=d)}{P(z_{i}<d)}\right)=\theta_{d}+\beta^{T}x_{i}$ (3)

where $\beta^{T}x_{i}$ is the propensity score for unit $u_{i}$and $d$ is the
treatment level. The parameter $\beta^{T}$ is estimated using maximum
likelihood to give the estimate $\hat{\beta}^{T}$. The propensity score can
now be incorporated into the edge weights of the graph. To further incentivise
the matching algorithm to pair units with the maximised difference in
treatments in the multiple treatment level case, their treatments are also
integrated into the weights. The distance metric is thus

$w_{ij}=\frac{(\hat{\beta}^{T}x_{i}-\hat{\beta}^{T}x_{j})^{2}+\epsilon}{(z_{i}-z_{j})^{2}}$
(4)

where $\epsilon$ is a very small positive real number put into place to
mitigate the edge case where the confounders of two units are the same which
invalidates the edge with a zero weight. With the graph and its edge weights
now defined, the optimal matching $M$ is found using a weight version of the
Edmond’s Blossom algorithm [25] which optimises Equation 2. Finally, the
average treatment effect is calculated as follows

$ATE=\frac{1}{M}\sum_{(u_{i},u_{j},w_{ij})\in
M}\frac{y_{i}-y_{j}}{z_{i}-z_{j}}$ (5)

##  IV Datasets

This section contextualises our problem by describing the observation data
collected and defining all the variables. The datasets were obtained from the
following two sources:

London Datastore222https://data.london.gov.uk/dataset/lsoa-atlas. The Greater
London Authority maintains an open data-sharing portal that contains census
data for each neighbourhood including statistics such as demographics,
employment, housing and transportation. The portal also provides neighbourhood
boundaries in the GeoJSON format.

Foursquare. The venue and tip data for the retail stores were obtained through
the LBSN platform Foursquare. Collected by Chen et al. [26], the dataset
consists of 19 million tips for 13 million venues from around the world. Each
venue includes information such as the venue ID, name, location and time of
account creation. Venues are also assigned categories such as “Restaurants”,
“Bank”, and “Pharmacy” and these categories have a hierarchical structure as
depicted in Figure 2. Tips are timestamped texts written by users for venues.

Figure 2: Hierarchy of Foursquare venue categories. Original hierarchy has
five levels and more categories in each level but only two are shown.

Since the geography of our problem is limited to the Greater London area, we
initially filtered venues using the country code field to shortlist all the
venues that belong to the United Kingdom. This significantly reduces the
number of venues from over 13 million to 272,079. These venues were assigned
to neighbourhoods based on whether any of the neighbourhood boundaries
enclosed their location (given by latitude and longitude). The unassigned
venues are consequently not in London and were discarded.

For reliability of results, the experiment is repeated thrice each year from
2011 to 2013. The treatment must temporally precede the outcome to measure the
causal effect. Hence, the treatment and outcome are measured at the beginning
and the end of the year respectively for each neighbourhood.

###  IV-A Treatment

The treatment variable must measure the presence of anchor stores which are
defined as large retail chains with a brand name [1]. Hence, a list of global
brand names was developed from all the values to the ‘brand’ key from
OpenStreetMaps 333https://taginfo.openstreetmap.org/keys/brand#values and the
names from the ‘chains’ field in Foursquare
444https://developer.foursquare.com/docs/build-with-foursquare/chains/. The
venues in Greater London are labelled as anchor stores if their name is
present in the brand name list and have the first level category of ‘Shop &
Service’. Only anchor stores that started before the year of the experiment
are counted. The treatment variable for each neighbourhood is the number of
anchor stores per restaurant (hence $Z\in\mathbb{R}^{+}$). To utilise the
propensity score described earlier, the treatment values are discretised into
three treatment levels. Since there is no significant literature on choosing
the number of treatment levels, we intuitively set it to three for two reasons
(1) The number of levels dictates the dimensionality of the propensity score
which must be set low; (2) Around 33% of the neighbourhoods do not have anchor
stores (i.e. treatment level is 0). Hence, we split the remainder of the
neighbourhoods into two groups yielding three different treatment levels so
that each level consists of roughly 33% of the neighbourhoods. Table I shows
the treatment value ranges for each treatment level.

###  IV-B Outcome

The outcome variable must measure the customer visit patterns in restaurants
in the Greater London area. Foursquare provides two types of customer
activities at the venues: check-ins and tips. Since a prior study [27]
indicated that 75% of the check-ins are fake, we use tips following [28]. Tips
provide a better proxy of customer visits because they must write a text on
the platform regarding the venue which serves as a verification of their
visit. The restaurants within each neighbourhood are determined using the
second-level categories containing the term ‘Restaurant’ for venues. Among
restaurants that started before the year of the experiment, the number of tips
posted on them during the year is counted. It is also ensured that the non-
anchor restaurants do not belong to a retail chain or brand name. Finally, the
outcome variable for each neighbourhood is the number of tips per restaurant.

TABLE I: Treatment level ranges in no. of anchor stores per restaurant in a neighbourhood year / treatment level | 0 | 1 | 2
---|---|---|---
2011 |  $\leq$ 0.000 |  $\leq$0.250 |  $>$0.250
2012 |  $\leq$0.000 |  $\leq$0.228 |  $>$0.228
2013 |  $\leq$0.000 |  $\leq$0.223 |  $>$0.223

Table II summarises the top 10 second-level categories and anchor brands.
Grocery stores are the most common type of anchor store which is consistent
with the literature [1, 8] which confirmed that large departmental stores
traditionally served as anchor stores. We observe that apparel stores,
electronics and pharmacies are also gaining attention as anchors.

TABLE II: Top 10 anchor store categories and brands by percentage of stores in Greater London till 2013 Category | % stores | Brand | % stores
---|---|---|---
Grocery Store | 23.67 | Tesco | 14.01
Clothing Store | 14.98 | Boots | 5.59
Pharmacy | 6.43 | Sainsbury’s | 5.54
Convenience Store | 5.33 | Co-op Food | 3.75
Bookstore | 4.27 | WHSmith | 3.29
Electronics Store | 3.08 | Iceland | 2.54
Hardware Store | 2.90 | Argos | 2.36
Shoe Store | 2.80 | Superdrug | 2.25
Office Supplies | 2.57 | Marks & Spencer | 2.08
Women’s Store | 2.48 | TK Maxx | 1.90

###  IV-C Confounding Variables

The selection of confounders is crucial to propensity score matching as
described earlier and requires domain knowledge of the retail industry. In our
context, confounding variables must be key factors for site selection that
impact the success of the retail stores and restaurants. To this end, we
consult several papers [29, 13, 30] and employ 15 confounders from the
following categories for each neighbourhood:

Population and Age. Retail stores not only seek populated neighbourhoods but
also target customers of a particular demographic. This is captured by three
variables measuring population for age groups of (i) 0-15, (ii) 16-64 (iii)
65+.

Education. Education plays an important role in determining the living
standard and availability of skilled labour pool which is measured using the
Average GCSE point scores.

Housing. The standard of living is also reflected through affordable housing
accounted for by the following variables (i) median house price, (ii) number
of dwellings (iii) the percentage of dwellings sold during the year.

Employment. The reason for including employment-related variables is two-fold.
Employment indicators not only suggest the availability of the labour pool but
also potential customers with disposal income. This is reflected by the
following three variables: (i) number of full-time (ii) number of part-time
employees (iii) their mean income.

Business Accessibility. This category measures the ability of customers to pay
a visit to the retail store or restaurant. Tayeen et al. [30] suggest that the
availability of transportation access and low crime rate are relevant
variables. These are measured by the average public transport accessibility
score (PTAL) and the total crime rate in the datastore.

Tourism. Several customers are attracted to a restaurant’s neighbourhood not
only because of the presence of anchor stores but also its tourism
significance. According to Tayeen et al., [30], the tourism significance of a
neighbourhood is captured based on the counts of four types of attractions:
(i) tourist locations (ii) rivers, lakes and reservoirs (iii) parks. These are
measured by counting venues in Foursquare with corresponding categories in
each neighbourhood of London.

##  V Results

This section presents the results of the experiments. We first match the
neighbourhoods based on their propensity scores. We then evaluate the
propensity score matching through three metrics namely: dose difference,
confounder balance and performance on a synthetic dataset [31]. We finally
compare the outcome of all matched pairs of neighbours and quantify the effect
of anchor stores on customer traffic in Greater London.

###  V-A Correlation Analysis

To obtain an initial insight, we plotted the number of stores per restaurant
(i.e. treatment) against the number of tips per restaurant (i.e. outcome)
which is shown in figure 3. The spearman correlation [32] is utilised to
measure the correlation between the two variables. A small but significant
positive correlation is observed between the treatment and outcome variables
across all three years. Although premature, the results suggest that a higher
number of anchor stores in a neighbourhood correlate with an increase in the
number of tips to restaurants. Yet, the results from the following sections
will establish a causal relationship between the two variables.

Figure 3: Correlation analysis between the treatment and outcome variables. TABLE III: Counts of dose levels between matched pairs of (low dose, high dose) units. year / dose pair | (0, 1) | (1, 2) | (0, 2)
---|---|---|---
2011 | 103 | 38 | 119
2012 | 93 | 66 | 121
2013 | 97 | 80 | 113

###  V-B Dose Difference

The dose difference exposes the efficacy of a specific aspect in the matching
method: maximising the difference in treatment levels between matched pairs.
This mitigates the influence of noise when calculating the average treatment
effect. The matrix plots, as shown in Table III, depict the distributions of
doses in matched pairs for the three years (i.e., 2011, 2012 and 2013). The
rows represent the treatment received by a unit at a higher level compared to
its match. The column represents the treatment level of the unit at a lower
level than its counterpart. The treatment difference is two levels for 45.8%,
43.2% and 39.0% of the pairs in 2011, 2012 and 2013 respectively. While the
remainder of the percentages reflects the pairs with a treatment level
difference of one.

TABLE IV: The average confounder values between high and low dose units among matched pairs in 2013. | high | low | t-test | p-value
---|---|---|---|---
% of dwellings sold | 3.29 | 3.22 | 0.79 | 0.43
Aged 0-15 | 2739.14 | 2647.41 | 1.32 | 0.19
Aged 16-64 | 9526.72 | 9216.72 | 1.81 | 0.07
Aged 65+ | 1521.03 | 1517.41 | 0.09 | 0.93
Average GCSE scores | 350.91 | 349.30 | 1.09 | 0.28
Average PTAL score | 3.93 | 3.78 | 1.33 | 0.18
Mean Income | 53676.62 | 52758.07 | 0.65 | 0.52
Median House Price | 391187.98 | 375510.38 | 0.78 | 0.43
No. of Full-time employees | 6179.66 | 5147.93 | 0.68 | 0.50
No. of Part-time employees | 2206.21 | 1871.72 | 1.06 | 0.29
No. of dwellings | 5639.12 | 5445.98 | 1.99 | 0.05
No. of lakes & rivers | 0.14 | 0.13 | 0.19 | 0.85
No. of parks | 1.19 | 1.16 | 0.23 | 0.82
No. of tourist locations | 1.31 | 1.67 | -0.69 | 0.49
Total crime rate | 98.77 | 90.65 | 0.94 | 0.35

###  V-C Confounder Balance

The matching algorithm is also evaluated for its effectiveness to pair units
with high confounder similarity. Table IV shows the average confounder values
for units with higher and lower treatment levels in each pair for the year
2013 (other years also hold similar results). The independence t-test was
performed to realise the significance of the confounder similarity. For every
confounder variable, the null hypothesis that its average value in high
treatment units of each pair is identical to the average value in low
treatment units of each pair cannot be rejected (i.e., the confounder values
are highly similar between paired units). Furthermore, the absolute value of
the test statistic is less than 1 for most confounders suggesting a
significant confounder similarity between the matched pairs of high and low
treatment levels.

Figure 4: Confounder distributions from observational data

###  V-D Treatment Effect

Now that it has been established that high treatment units and low treatment
units have been sampled from almost identical populations, the average
treatment effect is now computed. Table V shows the average treatment effects
calculated across all three years. The ATE, in this case, indicates the change
in the number of tips per anchor store across all restaurants in the store’s
neighbourhood. An ATE value of zero implies no effect. Each anchor store in
the neighbourhood caused an average increase between 0.05 and 0.15 tips to all
restaurants in the neighbourhood. On average, this amounts to a 14.2-26.5%
increase in tips per restaurant between the matched neighbourhoods per year.
Moreover, for all three years, ATE is positively skewed suggesting that anchor
stores in a neighbourhood cause customer traffic to non-anchor restaurants.

TABLE V: Distribution of treatment effects across three years. year | ATE | min | 25% | 75% | max
---|---|---|---|---|---
2011 | 0.13 | -5.10 | -0.50 | 0.75 | 9.53
2012 | 0.19 | -3.33 | -0.33 | 0.56 | 4.95
2013 | 0.16 | -5.60 | -0.38 | 0.58 | 35.48

###  V-E Synthetic Dataset

An effective way to evaluate the propensity score matching with known outcomes
is to develop a synthetic dataset, following several works in this area [22,
33, 34]. However, rather than developing an entirely random dataset, we
generate the synthetic data with supervision from observed data. In essence,
the influence of confounders on the treatment and the outcome as well as the
treatment’s influence on the outcome is modelled. Hence, the true causal
effect of the treatment is known beforehand. Then, the matching method’s
efficacy to recover the true effect can be evaluated.

To achieve this, we first generate $N$ units by generating confounders for
each unit. The confounders are generated by modelling each confounder’s
distribution and sampling from it $N$ times. Based on the distributions of the
confounders from the observed data as shown in Figure 4, it can be discerned
that truncated normal distribution is used to model the confounders except for
the number of full-time and part-time employees and the three tourism
variables which follows the exponential distribution. Maximum likelihood
estimation was utilised to determine the parameters of the distributions.

Then, two functions are modelled: $f_{Z}(x_{i})$ generates the treatment given
the confounders, while $f_{Y}(z_{i},x_{i})$ generates the outcome given the
confounders and treatment. Based on the distributions of treatment and
outcomes as depicted in figure 5, exponential distribution can model the
treatment as:

$f_{Z}(x)=P(Z=z|x_{i})=\begin{cases}\frac{1}{f(x_{i})}e^{\frac{z}{f(x_{i})}},&z\geq
0\\\ 0,&z<0\\\ \end{cases}$ (6)

$f(x_{i})=\frac{\sum_{p}x_{ip}}{\eta}$ (7)

The truncated normal distribution is used to determine the outcome for each
unit:

$f_{Y}(z_{i},x_{i})=Y_{i}\sim
Norm\left(\alpha\sum_{p}x_{ip}+\beta\sum_{p}x_{ip}\right)+\gamma z_{i}$ (8)

$\alpha$, $\beta$ and $\eta$ are tuned until the range of values matches those
of the observed data. $\gamma$ represents the value of the ground truth
average treatment effect. Figure 5 shows the treatment and outcome
distributions for the observational and synthetic datasets are highly similar
when these parameters are tuned.

The synthetic data was generated for four values of average treatment effect.
A total of 6250 units were generated for each treatment effect, while the
matching was run on 10 batches of 625 units. The treatment effects values
include large and small positive and negative values. Table VI summarises the
results. The estimated ATEs from synthetic data are highly similar to the
ground truth ATEs. While the treatment effects take a range of values across
all cases, at least 75% of effects are positive effects in most cases. This is
consistent with the results obtained from the real data. We also observe that
the matching performs well on larger ground truth ATE compared to the lower
values as the difference between true and estimated ATE is large for smaller
ATE values. This suggests that one must take a cautious stance when
interpreting the results with real data in Table 6.

Figure 5: Top: true treatment and outcome distributions, Bottom: generated treatment and outcome distributions TABLE VI: Treatment effects for synthetic data generated using various ground truth average treatment effects.  ATE | est. ATE | min | 25% | 75% | max
---|---|---|---|---|---
0.25 | 0.22 | -57.11 | -0.84 | 1.22 | 48.52
1.50 | 2.06 | -74.76 | 0.52 | 2.57 | 1430.55
4.00 | 4.14 | -41.69 | 3.03 | 5.12 | 170.27
10.00 | 9.80 | -71.70 | 8.92 | 11.02 | 60.46

Based on the results from correlation analysis and matching, we observe that
neighbourhood characteristics have a significant impact on the presence of
anchor stores and their subsequent causation of increased customer visitation
to non-anchor stores. The positive treatment effects observed across all three
years reinforce previous work in Economics on the anchor store effect. This
consistency in results with the literature and the performance across the
presented evaluation metrics reveal the viability of applying causal inference
to investigations delving into the retail industry.

##  VI Discussion and Future Work

The above findings lead us to discuss the implications and provide practical
application scenarios for researchers, small business owners, mall developers,
and mall managers. A few examples are as follows.

  1. ($\mathcal{I}_{1}$)

Non-anchor stores can utilise the causal analysis for site selection decisions
by targeting neighbourhoods with significant anchor store effects. For
instance, they can compare the customer traffic between current and new
locations based on the anchor store presence.

  2. ($\mathcal{I}_{2}$)

Similarly, The low dose units in matched pairs with a dose difference serve as
future location recommendations to anchor stores and mall developers as they
demonstrate similar neighbourhood conditions to the high dose counterparts
with lower competition.

  3. ($\mathcal{I}_{3}$)

Mall managers can derive the value of their malls from the anchor store
effects in their neighbourhoods and develop lucrative rental contracts for
their tenants.

Nonetheless, there are limitations in our work which must be mentioned along
with potential future works.

Unknown Confounding Variables. In our analysis, we took into account 15
confounding variables of 6 different aspects. Nevertheless, there may still
exist several variables with confounding effects for example various social
events that impact customer traffic to anchor and non-anchor stores alike. As
pointed out in [31], uncertainty pertains to whether all confounding variables
have been explored. A systematic approach to solving this challenge remains an
open problem.

Data Availability. The dataset studied is bottlenecked by time and space due
to the limited data availability. The census data below annual granularity and
for other cities has been challenging to find. Also, the dataset only covers
three years, limiting our analysis to the period from 2011 to 2013. The
geography and time in our study are discretised by neighbourhood and year
respectively, which might be too broad. It is worthwhile to employ causal
modelling like Granger Causality [35] to leverage spatiotemporal properties of
LBSNs and study this effect at a more fine-grained level [36, 37].

##  VII Conclusion

In this paper, we adopted a causal approach and quantitatively demonstrated
the anchor store effect. Contrary to previous work in this area that perform
studies on confidential datasets, we employed publicly available data from an
LBSN and a government’s data-sharing portal making our study more cost-
effective. Additionally, we employed a causal analysis framework based on
propensity score matching to mitigate the bias induced by confounding effects
between variables that impact customer traffic to anchor stores and non-anchor
stores. As a case study, we applied this framework to over 600 neighbourhoods
in the Greater London area and concluded that the anchor stores cause an
increase of 14.2-26.5% tips per year to non-anchor stores within their
neighbourhoods. Furthermore, we performed extensive evaluations that
demonstrated the framework’s high matching quality and resilience to noise
through experiments on synthetic datasets with known treatment effects.

Though we have realised a causal relation in its simplest form between two
variables (i.e. treatment and outcome), we hope that further investigations
will be undertaken to discern more complex causal mechanisms in Retail
Economics.

## Acknowledgments

This research has been supported in part by the 5GEAR project (Decision No.
318927) and the FIT project (Decision No. 325570) from the Academy of Finland.

## References

  * [1] H. Konishi and M. T. Sandfort, “Anchor stores,” _Urban Economics_ , 2003.
  * [2] M. A. Vitorino, “Empirical entry games with complementarities: An application to the shopping center industry,” _Marketing Research_ , vol. 49, 2012.
  * [3] M. J. Eppli and J. D. Benjamin, “The evolution of shopping center research: A review and analysis,” _Real Estate Research_ , vol. 9, 1994.
  * [4] S. D. Diana, D. C. José, and C. P. José, “The impact of anchor stores on the performance of shopping centres: the case of Sonae Sierra,” _International Retail & Distribution Management_, vol. 39, 2011.
  * [5] B. P. Pashigian and E. D. Gould, “Internalizing externalities: the pricing of space in shopping malls,” _The Law & Economics_, vol. 41, 1998.
  * [6] J. D. Benjamin, G. W. Boyle, and C. F. Sirmans, “Retail Leasing: The Determinants of Shopping Center Rents,” _AREUEA_ , vol. 18, 1990.
  * [7] C. Sirmans and K. Guidry, “The determinants of shopping center rents,” _Real Estate Research_ , vol. 8, 1993.
  * [8] A. Finn and J. J. Louviere, “Shopping center image, consideration, and choice: Anchor store contribution,” _Business Research_ , vol. 35, 1996.
  * [9] P. M. Anderson, “Assoc. of shopping-center anchors with performance of a nonanchor specialty chains stores,” _Retailing_ , vol. 61, 1985.
  * [10] S. M. Davidson and A. Rummel, “Retail changes associated with wal-mart’s entry into maine,” _International Retail & Distribution Management_, 2000\.
  * [11] G. M. Artz and K. E. Stone, “Analyzing the impact of wal-mart supercenters on local food store sales,” _American agricultural economics_ , vol. 88, 2006\.
  * [12] T. Hernandez, “The impact of big box internationalization on a national market: a case study of home depot inc. in canada,” _The International Review of Retail, Distribution and Consumer Research_ , vol. 13, 2003.
  * [13] S.-O. Daunfeldt, O. Mihaescu, H. Nilsson, and N. Rudholm, “Spillover effects when ikea enters: Do incumbent retailers win or lose?” _Papers in Regional Science_ , vol. 98, 2019.
  * [14] P. R. Rosenbaum and D. B. Rubin, “The central role of the propensity score in observational studies for causal effects,” _Biometrika_ , vol. 70, 1983\.
  * [15] D. B. Rubin, “Estimating causal effects of treatments in randomized and nonrandomized studies,” _Educational Psychology_ , vol. 66(5), 1974.
  * [16] T. C. Chalmers, H. Smith Jr, B. Blackburn, B. Silverman, B. Schroeder, D. Reitman, and A. Ambroz, “A method for assessing the quality of a randomized control trial,” _Controlled clinical trials_ , vol. 2, 1981.
  * [17] J. Pearl _et al._ , “Causal inference in statistics: An overview,” _Statistics surveys_ , vol. 3, 2009.
  * [18] D. L. Thistlethwaite and D. T. Campbell, “Regression-discontinuity analysis: An alternative to the ex post facto experiment.” _Educational psychology_ , vol. 51, 1960.
  * [19] D. Acemoglu, S. Johnson, and J. A. Robinson, “The colonial origins of comparative development: An empirical investigation,” _American economic review_ , vol. 91, 2001.
  * [20] P. R. Rosenbaum and D. B. Rubin, “The central role of the propensity score in observational studies for causal effects,” _Biometrika_ , vol. 70, 1983.
  * [21] B. Lu, E. Zanutto, R. Hornik, and P. R. Rosenbaum, “Matching with doses in an observational study of a media campaign against drug abuse,” _the American Stat. Assoc._ , vol. 96, 2001.
  * [22] S. Li, N. Vlassis, J. Kawale, and Y. Fu, “Matching via dimensionality reduction for estimation of treatment effects in digital marketing campaigns,” in _IJCAI_ , 2016.
  * [23] M. M. Joffe and P. R. Rosenbaum, “Invited Commentary: Propensity Scores,” _American Epidemiology_ , 1999.
  * [24] P. McCullagh, “Regression models for ordinal data,” _the Royal Stat. Society. Series B (Methodological)_ , 1980.
  * [25] Z. Galil, “Efficient algorithms for finding maximum matching in graphs,” _ACM Computing Surveys_ , 1986.
  * [26] Y. Chen, J. Hu, H. Zhao, Y. Xiao, and P. Hui, “Measurement and analysis of the swarm social network with tens of millions of nodes,” _IEEE Access_ , vol. 6, 2018.
  * [27] Z. Zhang, L. Zhou, X. Zhao, G. Wang, Y. Su, M. Metzger, H. Zheng, and B. Y. Zhao, “On the validity of geosocial mobility traces,” in _Proc. of the Twelfth ACM Workshop on Hot Topics in Networks_ , 2013.
  * [28] Y. D. Kwon, D. Chatzopoulos, E. ul Haq, R. C.-W. Wong, and P. Hui, “Geolifecycle: User engagement of geographical exploration and churn prediction in lbsns,” _IMWUT_ , 2019.
  * [29] E. M. Mazze, “Identifying the key factors in retail store location,” _Small Business Management (pre-1986)_ , vol. 10, p. 17, 1972.
  * [30] A. S. M. Tayeen, A. Mtibaa, and S. Misra, “Location, location, location! quantifying the true impact of location on business reviews using a yelp dataset,” in _ASONAM_ , 2019.
  * [31] A. Hasthanasombat and C. Mascolo, “Understanding the Effects of the Neighbourhood Built Environment on Public Health with Open Data,” in _WWW_ , 2019.
  * [32] E. C. Fieller, H. O. Hartley, and E. S. Pearson, “Tests for rank correlation coefficients. i,” _Biometrika_ , vol. 44, 1957.
  * [33] X. S. Gu and P. R. Rosenbaum, “Comparison of multivariate matching methods: Structures, distances, and algorithms,” _Computational and Graphical Statistics_ , 1993.
  * [34] F. Tsapeli, M. Musolesi, and P. Tino, “Non-parametric causality detection: An application to social media and financial data,” _Physica A: Stat. Mechanics and its Applications_ , vol. 483, 2017.
  * [35] C. W. J. Granger, “Investigating causal relations by econometric models and cross-spectral methods,” _Econometrica_ , vol. 37, 1969.
  * [36] A. K. Vallapuram, N. Nanda, Y. D. Kwon, and P. Hui, “Interpretable business survival prediction,” in _ASONAM_ , 2021.
  * [37] L.-y. Chen, Y. Chen, Y. D. Kwon, Y. Kang, and P. Hui, “Ian: Interpretable attention network for churn prediction in lbsns,” in _ASONAM_ , 2021.


Stochastic gravitational wave background
from early dark energy

Naoya Kitajima${}^{\,a,b}$ and Tomo Takahashi${}^{\,c}$

aFrontier Research Institute for Interdisciplinary Sciences, Tohoku
University, Sendai 980-8578, Japan
bDepartment of Physics, Tohoku University, Sendai 980-8578, Japan
cDepartment of Physics, Saga University, Saga 840-8502, Japan

We study the production of stochastic gravitational wave background from early
dark energy (EDE) model. It is caused by resonant amplification of scalar
field fluctuations, which easily takes place for typical EDE potential based
on the string axion or $\alpha$-attractor model. The resultant spectrum of
gravitational wave background is computed by performing 3D lattice
simulations. We show that, specifically in some class of generalized
$\alpha$-attractor EDE model, a significant amount of gravitational waves can
be produced via tachyonic instability with a peak around femto-Hz frequency
range. Models predicting such gravitational waves can be constrained by the
cosmic microwave background observations.

##  1 Introduction

Recent observational discrepancy of the Hubble constant $H_{0}$, known as the
Hubble tension, where the values of $H_{0}$ derived from direct and indirect
measurements are inconsistent at almost 5$\sigma$ level (see e.g., [1, 2] for
a review of the current status of the tension), motivates us to consider the
new physics beyond the $\Lambda$CDM model. One of such models suggested to
resolve the Hubble tension is an early dark energy (EDE) scenario [3] in which
a scalar field fractionally contributes to the total energy density at around
recombination epoch. As a potential for the early dark energy, the axion type
one is assumed in many works. Indeed axion is one of the representative
example of the new physics which is originally motivated by the particle
physics [4, 5, 6, 7] and is also predicted by string theory [8, 9, 10].
Therefore axion can be a plausible candidate for the EDE because it has
generally a large initial misalignment amplitude while the mass is extremely
small [11, 12]. An extremely small mass is required because the axion should
be frozen until the recombination epoch to be regarded as EDE and a large
initial amplitude is also needed because the axion should give a sizable
contribution to the total energy density in order to affect the expansion of
the Universe, which can reduce the sound horizon to increase the value of
$H_{0}$ when fitted to cosmic microwave background (CMB) data such as Planck
[13].

Another motivation for the axion EDE is the cosmic birefringence, recently
reported in [14] using the Planck data. In general, the background axion EDE
field can induce the rotation of the linear polarization angle of the CMB
photon and hence it can explain the reported isotropic cosmic birefringence
[15, 16, 17, 18, 19].

In this paper, we discuss the consequences of the axion EDE other than the
Hubble tension and the cosmic birefringence. In particular, we show that the
stochastic gravitational wave (GW) background is predicted from the axion EDE
scenario through the 2nd order scalar field fluctuations. Actually the scalar
field fluctuations can be resonantly amplified through the self-interaction or
interaction with other field as in the case of the preheating after inflation
[20, 21]. For the case of the preheating, a significant amount of GW can be
produced with the spectrum peaked typically at GHz range [22, 23, 24, 25, 26,
27, 28]. By the same mechanism, string axion can also induce stochastic GWs in
lower frequency bands depending on the axion mass [29, 30, 31, 32] and axion
dark matter scenario typically predicts GWs with nHz frequency range [33, 34,
35, 36, 37, 38]. Here in this paper, we argue that the GWs are inevitably
generated in the axion EDE scenario for some particular parameter choice [39],
which can be regarded as a unique cosmological signature of the axion
EDE#1#1#1 Actually GW production from the axion EDE model is studied in [40]
where the explosive gauge-field production is assumed, which is a different
scenario from the one we consider in this paper. . To this end, we assume two
types of potential for the EDE, which has been motivated by string theories
and show that the GW spectrum is peaked around the frequency of $f\sim
10^{-15}~{}{\rm Hz}$ which is a new source of GWs at the frequency region.

The structure of this paper is as follows. In the next section, we describe
the axionic early dark energy model we consider in this paper and discuss the
evolution of the scalar field and the spectrum of density fluctuations
obtained from lattice simulation. Then in Section 3, we discuss the production
of GWs in the scenario. The GW spectrum generated in our model is presented.
The final section is devoted to conclusions and discussions.

##  2 Early dark energy and parametric resonance

Figure 1:  The shape of the potential for the axion EDE type with $n=2$
(dotted), the $\alpha$-attractor-type with $(p,q,r,\beta)=(2,0,0,0)$ (dashed)
and the generalized $\alpha$-attractor-type with $(p,q,r,\beta)=(2,2,1,10)$
(solid). For the $\alpha$-attractor cases, $f$ corresponds to
$f=\sqrt{6\alpha}M_{P}$.

Let us consider the model with the following axion-like EDE potential,

$\displaystyle
V(\phi)=m^{2}f^{2}\left[1-\cos\bigg{(}\frac{\phi}{f}\bigg{)}\right]^{n},$ (1)

where $n$ is an integer, $f$ is the decay constant and $m$ is the energy scale
which determines the onset of the field oscillation. We also consider a
generalized $\alpha$-attractor type potential [41, 42] as follows:

$\displaystyle
V(\phi)=\frac{V_{0}(1+\beta)^{2r}\tanh^{2p}(\phi/\sqrt{6\alpha}M_{P})}{[1+\beta\tanh^{q}(\phi/\sqrt{6\alpha}M_{P})]^{2r}}\,,$
(2)

where $p$, $q$ and $r$ are (positive) integers and $\beta$ is a numerical
constant. $V_{0}$ corresponds to the potential value at $\phi\to\infty$.
$\alpha$ is a parameter which controls the energy scale of the scalar field.
Both of the above potentials can be approximated by
$m^{2}f^{2}(\phi/f)^{2n}/2^{n}$ (axionlike), and
$V_{0}(1+\beta)^{2r}(\phi/\sqrt{6\alpha}M_{P})^{2p}$ ($\alpha$-attractor) for
small $\phi$ limit. The shape of the potential is illustrated in Figure 1.
Assuming the flat Friedmann‐Lemaître‐Robertson-Walker (FLRW) background
metric, the field equation is given by

$\displaystyle\ddot{\phi}+3H\dot{\phi}-\frac{\nabla^{2}\phi}{a^{2}}+\frac{\partial
V}{\partial\phi}=0,$ (3)

where the overdot represents the derivative with respect to the cosmic time,
$a$ is the scale factor and $H=\dot{a}/a$ is the Hubble parameter. The Hubble
parameter is determined by the background radiation/matter mixed fluid. The
scalar field is initially placed far from the minimum (i.e. misaligned) and
almost frozen due to the Hubble friction term. During this period, the energy
density of the scalar field is kept constant and thus it behaves like a dark
energy with equation of state $w\simeq-1$. The scalar field starts to roll
down the potential when the Hubble friction becomes inefficient against the
potential gradient and then it oscillates around the minimum of the potential.
In this phase, it is no longer regarded as dark energy, and its energy density
dilutes faster than that of matter component in order to be consistent with
the observed $\Lambda$CDM cosmology when the EDE field starts to oscillate at
early times. Therefore, the $n=1$ case is in general not compatible with the
early dark energy model especially in the context of the Hubble tension#2#2#2
Considering the interaction between the axion and dark photon, the axion
abundance can be significantly suppressed due to the tachyonic production of
dark photons. The suppression factor can be typically as small as $10^{-2}$ –
$10^{-3}$ [43, 44, 35]. It may open up a possibility of the axion EDE with
$n=1$. .

The oscillation of the scalar field is highly anharmonic in the case with
$n\geq 2$. In other words the self-interaction always exists. In this case,
the nonzero (inhomogeneous) mode can be excited as a consequence of tachyonic
instability or parametric resonance instability. It shows exponential
amplifications of the field fluctuation and the inhomogeneous mode easily
dominates over the homogeneous coherent oscillation mode. In particular, the
case with $n=2$ is a special one because it approaches conformally invariant
model ($V=(\lambda/4)\phi^{4}$) in the small field limit. In this case, the
resonant amplification persistently occurs until the backreaction turns on
[45]. In the subsequent sections, we show the dynamics of this amplification
process by numerical lattice simulations.

In what follows, we focus on the following specific cases:

  * •

(a) Axion EDE-type potential with $n=2$,

  * •

(b) Simple $\alpha$-attractor-type potential with $(p,q,r,\beta)=(2,0,0,0)$,

  * •

(c) Generalized $\alpha$-attractor-type potential with
$(p,q,r,\beta)=(2,2,1,10)$,

  * •

(d) Generalized $\alpha$-attractor-type potential with
$(p,q,r,\beta)=(3,2,1.5,10)$.

###  2.1 Lattice simulation

Because the scalar field fluctuation (i.e. the nonzero mode) grows
exponentially due to the tachyonic/parametric resonance instability, it soon
becomes comparable to the homogeneous mode and the growth is saturated at that
time. Then, the produced nonzero mode start to backreact to the homogeneous
mode and the dynamics enters the nonlinear regime. Thus, in order to correctly
follow the full dynamics of the system, we need to directly solve the
evolution equation (3) with discretized lattice space where the spatial
derivative is replaced with the finite difference. The setup of our simulation
for each case is summarized in Table 1. The initial field value for the
homogeneous mode is $\phi_{i}/f=3$ in all cases ($f=\sqrt{6\alpha}M_{P}$ for
the $\alpha$-attractor cases) and we put random Gaussian field on each grid
point with the standard deviation $\delta\phi=10^{-5}$ as the initial field
fluctuations. The choices of other parameters in the potential have also been
listed in Table 1. With the parameter values assumed here, the EDE can give a
sizable energy density fraction at around the matter-radiation equality.
Indeed a lot of analysis using observational data such as CMB have been
performed to find the favored values for the EDE parameters (see, e.g., [46]
and the references therein). We may take the best-fit values from those
analysis to choose the values of the EDE parameters, however, the preferred
ones can depend on the data set adopted and the setup in the analysis. We
should also mention that the nonlinear effect discussed below could affect the
evolution of EDE to change the fit to the data. In the light of these
considerations, we just assume the values listed in Table 1 as examples. We
show the time evolution of the EDE field from our simulations in Fig. 2. The
red line shows the time evolution of the scalar field, which exhibits the
coherent oscillation and the blue line represents the root-mean-square of the
density fluctuation of EDE, $\delta\rho_{\rm EDE}/\rho_{\rm EDE}$. The green
line corresponds to the energy fraction of EDE to the background energy
density, $\rho_{\rm EDE}/(\rho_{r}+\rho_{m})$ with $\rho_{r}$ and $\rho_{m}$
being that of radiation and matter, respectively. The dashed green line
corresponds to the energy fraction of EDE for the homogeneously oscillating
case without the nonlinear effect (backreaction). Interestingly, by comparing
the energy density with and without the nonlinear effect, one can clearly see
that the backreaction modifies the evolution of the EDE energy density, which
may affect the constraint on EDE from cosmological data such as CMB. Although
this would be an important issue, we here focus on the production of GWs
through the 2nd order scalar field fluctuations.

(a) Axion EDE

(b) $\alpha$-attractor EDE

(c) Generalized $\alpha$-attractor ($p=2$)

(d) Generalized $\alpha$-attractor ($p=3$)

Figure 2:  Time evolution of the zero mode oscillation (red), energy fraction of EDE (green) and root-mean-square of density fluctuation (blue) of EDE. The horizontal axis represents the redshift. The dashed green line corresponds to the energy fraction of EDE including homogeneous component (without backreaction).  | grid points | initial box size |  $z$-range | $f$ | $m$ | $\alpha$ | $V_{0}$
---|---|---|---|---|---|---|---
(a) | $512^{3}$ | $0.8\pi H_{\rm eq}^{-1}$ |  $10^{4}$ – $1$ | $0.1M_{P}$ | $5H_{\rm eq}$ | - | -
(b) | $512^{3}$ | $0.4\pi H_{\rm eq}^{-1}$ |  $10^{4}$ – $10$ | - | - | $10^{-4}$ | $0.17H_{\rm eq}^{2}M_{P}^{2}$
(c) | $512^{3}$ | $0.04\pi H_{\rm eq}^{-1}$ |  $10^{4}$ – $500$ | - | - | $10^{-5}$ | $0.3H_{\rm eq}^{2}M_{P}^{2}$
(d) | $512^{3}$ | $0.04\pi H_{\rm eq}^{-1}$ |  $10^{4}$ – $500$ | - | - | $10^{-5}$ | $0.3H_{\rm eq}^{2}M_{P}^{2}$
Table 1: Setup of our lattice simulations. Here $H_{\rm eq}$ is the Hubble
parameter at the time of matter-radiation equality.

For the axion EDE case (case (a)), the axion starts to oscillate at the
redshift $z_{\rm osc}\sim 4000$ (which corresponds to $z_{c}$, the redshift at
which the density fraction of EDE is maximized), but the growth of scalar
field fluctuations starts well after the onset of the oscillation ($z\sim 60$)
as also shown by the linear analysis in [47]. The growth stops when the
density contrast becomes ${\cal O}(1)$ at $z\sim 10$. For a simple
$\alpha$-attractor model (case (b)), the growth starts earlier (at $z\sim
600$) and saturates at $z\sim 100$. Note that in the above two cases (a) and
(b), the resonant amplification starts after the last scattering surface of
CMB ($z\sim 1000$) and hence it has no chance to affect the CMB observations
such as the B-mode signal due to the gravitational waves.

On the other hand, one can see that efficient amplification happens in the
cases (c) and (d). It is due to the tachyonic instability soon after the
commencement of the oscillation#3#3#3 The linear analysis of the field
fluctuation from the tachyonic instability with plateau type potential (like a
generalized $\alpha$-attractor model in our case) has been studied in [29, 30,
48]. . In this case, the growth rate is large (compared with the above two
cases) and the fluctuation becomes ${\cal O}(1)$ after several oscillations,
well before the last scattering surface. It should be noted that the efficient
amplification occurs even for $p=3$ case (case (d)).

(a) Axion EDE

(b) $\alpha$-attractor EDE

(c) Generalized $\alpha$-attractor ($p=2$)

(d) Generalized $\alpha$-attractor ($p=3$)

Figure 3:  Time evolution of the spectrum of the density fluctuation. The time
evolves from bottom to top and the thick black line corresponds to the
spectrum at the redshift $z=1$ (top-left), $10$ (top-right), and $100$ (bottom
panels). The horizontal axis is the wavenumber normalized by the horizon scale
at the matter-radiation equality, $k_{\rm eq}=a_{\rm eq}H_{\rm eq}$.

Fig. 3 shows the evolution of the spectrum of energy density of the EDE as a
function of wavenumber, which is normalized by $k_{\rm eq}=a_{\rm eq}H_{\rm
eq}$, corresponding to the horizon scale at the matter-radiation equality. The
time evolves from bottom to top and the thick black line corresponds to the
spectrum at redshift $z=10$ (top-left), $z=100$ (top-right) and $z=1000$
(bottom panels). In the top two figures 3(a) and 3(b) (axion EDE and simple
$\alpha$-attractor EDE models), one can see that the growth occurs in multiple
narrow bands and the spectrum has sharp peaks in the early stages. It is a
typical behaviour of the narrow resonance (see [30] for some detailed
discussion in a similar setup). After the growth is saturated, the spectrum is
smoothed through the momentum transfer due to the rescattering processes. In
the bottom two panels 3(c) and 3(d) (generalized $\alpha$-attractor models),
the growth is due to the tachyonic instability and thus the low-$k$ modes are
broadly enhanced in a short period of time. See [30] for some detailed
discussion regarding the growth due to the narrow resonance and the tachyonic
instability in a similar setup. Then, the nonlinear processes smooth out the
spectrum. Note that the spectrum has jagged shape in the top left figures
(case (a)). It is because the momentum transfer is not completed within the
simulation time. This is in contrast to other cases where the spectral shape
is smoothed due to the completion of the momentum transfer. Thus, the
nonlinear process can be imprinted in the temporal shape of the spectrum.

##  3 Gravitational wave production

(a) Axion EDE

(b) $\alpha$-attractor EDE

(c) Generalized $\alpha$-attractor ($p=2$)

(d) Generalized $\alpha$-attractor ($p=3$)

Figure 4:  Time evolution of the spectrum of gravitational waves. The time
evolves from bottom to top and the thick black line corresponds to the final
time, $z=1$ (top-left), $10$ (top-right), $100$ (bottom panels). The
horizontal axis is the wave number normalized by the horizon scale at the
matter-radiation equality. We have taken the initial value $\phi_{i}/f=3$.
Figure 5:  Spectrum of the GW density parameter at $z=1000$ for cases (c) and
(d) (solid lines), $100$ for case (b) (dashed line), and $10$ for case (a)
(dotted line). We have taken the initial value
$\phi_{i}/f=\phi_{i}/\sqrt{6\alpha}M_{P}=3$. Thin lines show the results with
larger box size ($0.08\pi H_{\rm eq}^{-1}$ with case (c) and (d)) and the gray
shaded region represents the current constraint from the CMB observation [49].

As shown in the previous section, fluctuations of the scalar EDE field can be
exponentially amplified as a consequence of the tachyonic and/or the
parametric resonance instabilities. Such a violent process leads to the
gravitational wave emission. The gravitational wave is defined by the tensor
metric perturbation $h_{ij}$ on the background FLRW metric which is given by

$\displaystyle ds^{2}=dt^{2}-a^{2}(t)(\delta_{ij}+h_{ij})dx^{i}dx^{j}.$ (4)

Here and in what follows, the sum is taken over repeated indices unless
otherwise stated. The evolution of gravitational waves is derived from the
linearlized Einstein equation:

$\displaystyle\ddot{h}_{ij}+3H\dot{h}_{ij}-\frac{\nabla^{2}h_{ij}}{a^{2}}=-16\pi
G\Lambda^{kl}_{ij}\Pi_{kl}$ (5)

where $\Lambda^{kl}_{ij}$ is the transverse-traceless projection tensor,
$\Pi_{ij}=-\partial_{i}\phi\partial_{j}\phi/a^{2}$ is the anisotropic stress
originated from the scalar field and $G$ is the Newton’s constant. The energy
density of stochastic GW is given by

$\displaystyle\rho_{\rm GW}(t)=\frac{1}{32\pi
G}\langle\dot{h}_{ij}(\bm{x})\dot{h}_{ij}(\bm{x})\rangle,$ (6)

where the angular bracket represents the ensemble average. Practically, this
average can be replaced with the spatial average in the lattice simulation.
The spectrum of GW density parameter is defined by

$\displaystyle\Omega_{\rm GW}(f)=\frac{1}{\rho_{\rm cr}}\frac{d\rho_{\rm
GW}}{d\ln f},$ (7)

where $\rho_{\rm cr}=3H^{2}M_{P}^{2}$ is the critical density and $f$ is the
frequency of GW.

We have solved the evolution equation for GWs (5) together with the scalar
field evolution (3) with the numerical lattice simulation following [26]. Here
we are interested in the induced GWs, and thus we neglect the primordial GW
background from the inflationary origin and set $h_{ij}=0$ as the initial
condition for our simulation.

The time evolution of the GW density parameter is shown in Fig. 4 for each
case. In the axion EDE case (figure 4(a)), the resonant amplification is not
efficient and thus the resultant GW amplitude (at $z=1$ shown by thick black
line) is smaller than that in other cases. The final GW amplitude is much
larger in a simple $\alpha$-attractor model (figure 4(b)), but still most GW
emission occurs only after the recombination epoch. Note that, in these cases,
the multiple-peak structure appears as an imprint of the narrow resonance band
structure. However, the emitted GWs cannot be probed by CMB in those cases. In
the cases (c) and (d) (figure 4(c) and 4(d)), low-$k$ modes are efficiently
excited due to the tachyonic instability, which induces large scalar
fluctuations. The amplitude of the emitted GW is many orders of magnitude
larger than the case with top panels#4#4#4 The density parameter of GWs looks
larger than 1 in this figure, but notice that the vertical axis is normalized
by $(f/M_{P})^{4}\sim 10^{-4}$, and hence $\Omega_{\rm GW}$ appears to be big,
however $\Omega_{\rm GW}$ does not exceed unity. . In this case, the efficient
GW emission occurs before the last scattering surface and it can affect the
CMB observation.

Figure 5 shows the resultant GW spectrum at $z=1000$ (solid line) for the
cases (c) and (d), 100 (dashed line) for the case (b), and 10 (dotted line)
for the case (a)#5#5#5At that time, the GW production is not yet completed and
thus the multiple-peak structure does not yet appear in the GW spectrum..
Although, due to the limitation of the lattice simulation, we have not
evaluated the GW spectrum at lower frequency, especially for the cases (c) and
(d), however if one simply extrapolates the solid lines to lower frequency,
e.g. $10^{-16}$ Hz, the CMB constraints may easily rule out the model [49].

##  4 Summary and discussion

In this paper, we have studied the GW emission from the EDE field fluctuations
amplified by the tachyonic/parametric resonance instabilities due to the self-
interactions. We have shown by performing the numerical lattice simulations
that the tachyonic instability efficiently enhances the field fluctuations and
the GWs can be produced before the recombination epoch while the narrow
resonance is less efficient and the GW emission is delayed. Note that the
tachyonic growth occurs near the minimum of the potential regardless of the
power index of the potential for the generalized $\alpha$-attractor type.

As discussed in Section 2, the nonlinear evolution of EDE may have a great
impact not only on the production of GWs, but also on the prospect of EDE as a
possible solution of the Hubble tension. In many works which study the fit of
the EDE to CMB and some other data, the evolutions of the background and
linear perturbation have been implemented in the analysis. However, we have
shown that the scalar field fluctuation can dominate over the coherent
oscillation. In such a case, the equation of state of the EDE can be modified
as in the preheating case [50]. Indeed, as shown in Fig. 2, the nonlinear
effect significantly reduces the energy density of the EDE compared with the
homogeneous case. As a consequence, it possibly changes the fit to CMB and
other cosmological data. Thus, one should correctly take into account the
temporal change of the equation of state of the EDE in the analysis to
investigate its constraints in some cases#6#6#6 The prediction for the cosmic
birefringence can also be affected. In particular, since the axion EDE field
becomes highly inhomogeneous, the anisotropic birefringence is expected rather
than the isotropic one and it can also be constrained by the CMB data [51,
52]. .

In this paper, we simply assume that the EDE does not have interaction with
other fields. As an extension of this work, one can consider a multi-field
model allowing the interaction between the EDE and other spectator fields (as
in the new EDE model [53, 54, 55]). In such a case, the resonant production of
spectator fields can be efficient and it can enhance the GW emission
accordingly. This issue would also be worth investigating, which is left for
future work.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant Number 19K03874 (TT), 20H01894
(NK), 20H05851 (NK), 21H01078 (NK), 21KK0050 (NK) and MEXT KAKENHI 23H04515
(TT).

## References

  * [1] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, “In the realm of the Hubble tension¥textemdasha review of solutions,” Class. Quant. Grav. 38 no. 15, (2021) 153001, arXiv:2103.01183 [astro-ph.CO].
  * [2] L. Perivolaropoulos and F. Skara, “Challenges for $\textyen Lambda$CDM: An update,” New Astron. Rev. 95 (2022) , arXiv:2105.05208 [astro-ph.CO].
  * [3] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, “Early Dark Energy Can Resolve The Hubble Tension,” Phys. Rev. Lett. 122 no. 22, (2019) 221301, arXiv:1811.04083 [astro-ph.CO].
  * [4] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [5] R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D16 (1977) 1791–1797.
  * [6] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [7] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [8] P. Svrcek and E. Witten, “Axions In String Theory,” JHEP 06 (2006) 051, arXiv:hep-th/0605206.
  * [9] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * [10] M. Cicoli, M. Goodsell, and A. Ringwald, “The type IIB string axiverse and its low-energy phenomenology,” JHEP 10 (2012) 146, arXiv:1206.0819 [hep-th].
  * [11] M. Kamionkowski, J. Pradler, and D. G. E. Walker, “Dark energy from the string axiverse,” Phys. Rev. Lett. 113 no. 25, (2014) 251302, arXiv:1409.0549 [hep-ph].
  * [12] T. Karwal and M. Kamionkowski, “Dark energy at early times, the Hubble parameter, and the string axiverse,” Phys. Rev. D 94 no. 10, (2016) 103523, arXiv:1608.01309 [astro-ph.CO].
  * [13] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. I. Overview and the cosmological legacy of Planck,” Astron. Astrophys. 641 (2020) A1, arXiv:1807.06205 [astro-ph.CO].
  * [14] Y. Minami and E. Komatsu, “New Extraction of the Cosmic Birefringence from the Planck 2018 Polarization Data,” Phys. Rev. Lett. 125 no. 22, (2020) 221301, arXiv:2011.11254 [astro-ph.CO].
  * [15] T. Fujita, K. Murai, H. Nakatsuka, and S. Tsujikawa, “Detection of isotropic cosmic birefringence and its implications for axionlike particles including dark energy,” Phys. Rev. D 103 no. 4, (2021) 043509, arXiv:2011.11894 [astro-ph.CO].
  * [16] G. Choi, W. Lin, L. Visinelli, and T. T. Yanagida, “Cosmic birefringence and electroweak axion dark energy,” Phys. Rev. D 104 no. 10, (2021) L101302, arXiv:2106.12602 [hep-ph].
  * [17] H. Nakatsuka, T. Namikawa, and E. Komatsu, “Is cosmic birefringence due to dark energy or dark matter? A tomographic approach,” Phys. Rev. D 105 no. 12, (2022) 123509, arXiv:2203.08560 [astro-ph.CO].
  * [18] K. Murai, F. Naokawa, T. Namikawa, and E. Komatsu, “Isotropic cosmic birefringence from early dark energy,” Phys. Rev. D 107 no. 4, (2023) L041302, arXiv:2209.07804 [astro-ph.CO].
  * [19] J. R. Eskilt, L. Herold, E. Komatsu, K. Murai, T. Namikawa, and F. Naokawa, “Constraint on Early Dark Energy from Isotropic Cosmic Birefringence,” arXiv:2303.15369 [astro-ph.CO].
  * [20] L. Kofman, A. D. Linde, and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett. 73 (1994) 3195–3198, arXiv:hep-th/9405187.
  * [21] L. Kofman, A. D. Linde, and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56 (1997) 3258–3295, arXiv:hep-ph/9704452.
  * [22] S. Y. Khlebnikov and I. I. Tkachev, “Relic gravitational waves produced after preheating,” Phys. Rev. D 56 (1997) 653–660, arXiv:hep-ph/9701423.
  * [23] R. Easther and E. A. Lim, “Stochastic gravitational wave production after inflation,” JCAP 04 (2006) 010, arXiv:astro-ph/0601617.
  * [24] R. Easther, J. T. Giblin, Jr., and E. A. Lim, “Gravitational Wave Production At The End Of Inflation,” Phys. Rev. Lett. 99 (2007) 221301, arXiv:astro-ph/0612294.
  * [25] J. Garcia-Bellido and D. G. Figueroa, “A stochastic background of gravitational waves from hybrid preheating,” Phys. Rev. Lett. 98 (2007) 061302, arXiv:astro-ph/0701014.
  * [26] J. Garcia-Bellido, D. G. Figueroa, and A. Sastre, “A Gravitational Wave Background from Reheating after Hybrid Inflation,” Phys. Rev. D 77 (2008) 043517, arXiv:0707.0839 [hep-ph].
  * [27] D. G. Figueroa and F. Torrenti, “Gravitational wave production from preheating: parameter dependence,” JCAP 10 (2017) 057, arXiv:1707.04533 [astro-ph.CO].
  * [28] P. Adshead, J. T. Giblin, and Z. J. Weiner, “Gravitational waves from gauge preheating,” Phys. Rev. D 98 no. 4, (2018) 043525, arXiv:1805.04550 [astro-ph.CO].
  * [29] J. Soda and Y. Urakawa, “Cosmological imprints of string axions in plateau,” Eur. Phys. J. C 78 no. 9, (2018) 779, arXiv:1710.00305 [astro-ph.CO].
  * [30] N. Kitajima, J. Soda, and Y. Urakawa, “Gravitational wave forest from string axiverse,” JCAP 10 (2018) 008, arXiv:1807.07037 [astro-ph.CO].
  * [31] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A. Stefanek, “Audible Axions,” JHEP 01 (2019) 053, arXiv:1811.01950 [hep-ph].
  * [32] C. S. Machado, W. Ratzinger, P. Schwaller, and B. A. Stefanek, “Gravitational wave probes of axionlike particles,” Phys. Rev. D 102 no. 7, (2020) 075033, arXiv:1912.01007 [hep-ph].
  * [33] W. Ratzinger and P. Schwaller, “Whispers from the dark side: Confronting light new physics with NANOGrav data,” SciPost Phys. 10 no. 2, (2021) 047, arXiv:2009.11875 [astro-ph.CO].
  * [34] R. Namba and M. Suzuki, “Implications of Gravitational-wave Production from Dark Photon Resonance to Pulsar-timing Observations and Effective Number of Relativistic Species,” Phys. Rev. D 102 (2020) 123527, arXiv:2009.13909 [astro-ph.CO].
  * [35] N. Kitajima, J. Soda, and Y. Urakawa, “Nano-Hz Gravitational-Wave Signature from Axion Dark Matter,” Phys. Rev. Lett. 126 no. 12, (2021) 121301, arXiv:2010.10990 [astro-ph.CO].
  * [36] W. Ratzinger, P. Schwaller, and B. A. Stefanek, “Gravitational Waves from an Axion-Dark Photon System: A Lattice Study,” SciPost Phys. 11 (2021) 001, arXiv:2012.11584 [astro-ph.CO].
  * [37] R. T. Co, K. Harigaya, and A. Pierce, “Gravitational waves and dark photon dark matter from axion rotations,” JHEP 12 (2021) 099, arXiv:2104.02077 [hep-ph].
  * [38] E. Madge, W. Ratzinger, D. Schmitt, and P. Schwaller, “Audible axions with a booster: Stochastic gravitational waves from rotating ALPs,” SciPost Phys. 12 no. 5, (2022) 171, arXiv:2111.12730 [hep-ph].
  * [39] M. C. Johnson and M. Kamionkowski, “Dynamical and Gravitational Instability of Oscillating-Field Dark Energy and Dark Matter,” Phys. Rev. D 78 (2008) 063010, arXiv:0805.1748 [astro-ph].
  * [40] Z. J. Weiner, P. Adshead, and J. T. Giblin, “Constraining early dark energy with gravitational waves before recombination,” Phys. Rev. D 103 no. 2, (2021) L021301, arXiv:2008.01732 [astro-ph.CO].
  * [41] F. X. Linares Cedeño, A. Montiel, J. C. Hidalgo, and G. Germán, “Bayesian evidence for $\alpha$-attractor dark energy models,” JCAP 08 (2019) 002, arXiv:1905.00834 [gr-qc].
  * [42] M. Braglia, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, and K. Koyama, “Unified framework for early dark energy from $\alpha$-attractors,” Phys. Rev. D 102 no. 8, (2020) 083513, arXiv:2005.14053 [astro-ph.CO].
  * [43] N. Kitajima, T. Sekiguchi, and F. Takahashi, “Cosmological abundance of the QCD axion coupled to hidden photons,” Phys. Lett. B 781 (2018) 684–687, arXiv:1711.06590 [hep-ph].
  * [44] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, “Relic Abundance of Dark Photon Dark Matter,” Phys. Lett. B 801 (2020) 135136, arXiv:1810.07188 [hep-ph].
  * [45] P. B. Greene, L. Kofman, A. D. Linde, and A. A. Starobinsky, “Structure of resonance in preheating after inflation,” Phys. Rev. D 56 (1997) 6175–6192, arXiv:hep-ph/9705347.
  * [46] V. Poulin, T. L. Smith, and T. Karwal, “The Ups and Downs of Early Dark Energy solutions to the Hubble tension: a review of models, hints and constraints circa 2023,” arXiv:2302.09032 [astro-ph.CO].
  * [47] T. L. Smith, V. Poulin, and M. A. Amin, “Oscillating scalar fields and the Hubble tension: a resolution with novel signatures,” Phys. Rev. D 101 no. 6, (2020) 063523, arXiv:1908.06995 [astro-ph.CO].
  * [48] H. Fukunaga, N. Kitajima, and Y. Urakawa, “Efficient self-resonance instability from axions,” JCAP 06 (2019) 055, arXiv:1903.02119 [astro-ph.CO].
  * [49] T. J. Clarke, E. J. Copeland, and A. Moss, “Constraints on primordial gravitational waves from the Cosmic Microwave Background,” JCAP 10 (2020) 002, arXiv:2004.11396 [astro-ph.CO].
  * [50] D. I. Podolsky, G. N. Felder, L. Kofman, and M. Peloso, “Equation of state and beginning of thermalization after preheating,” Phys. Rev. D 73 (2006) 023501, arXiv:hep-ph/0507096.
  * [51] A. Gruppuso, D. Molinari, P. Natoli, and L. Pagano, “Planck 2018 constraints on anisotropic birefringence and its cross-correlation with CMB anisotropy,” JCAP 11 (2020) 066, arXiv:2008.10334 [astro-ph.CO].
  * [52] M. Bortolami, M. Billi, A. Gruppuso, P. Natoli, and L. Pagano, “Planck constraints on cross-correlations between anisotropic cosmic birefringence and CMB polarization,” JCAP 09 (2022) 075, arXiv:2206.01635 [astro-ph.CO].
  * [53] F. Niedermann and M. S. Sloth, “New early dark energy,” Phys. Rev. D 103 no. 4, (2021) L041303, arXiv:1910.10739 [astro-ph.CO].
  * [54] F. Niedermann and M. S. Sloth, “Resolving the Hubble tension with new early dark energy,” Phys. Rev. D 102 no. 6, (2020) 063527, arXiv:2006.06686 [astro-ph.CO].
  * [55] J. S. Cruz, S. Hannestad, E. B. Holm, F. Niedermann, M. S. Sloth, and T. Tram, “Profiling Cold New Early Dark Energy,” arXiv:2302.07934 [astro-ph.CO].