# Wiki-based Communities of Interest: Demographics and Outliers

Hiba Arnaout1 Simon Razniewski1, Jeff Z. Pan2

###### Abstract

In this paper, we release data about demographic information and outliers of
communities of interest. Identified from Wiki-based sources, mainly Wikidata,
the data covers 7.5k communities, such as members of the White House
Coronavirus Task Force, and 345k subjects, e.g., Deborah Birx. We describe the
statistical inference methodology adopted to mine such data. We release
subject-centric and group-centric datasets in JSON format, as well as a
browsing interface. Finally, we forsee three areas this research can have an
impact on: in social sciences research, it provides a resource for demographic
analyses; in web-scale collaborative encyclopedias, it serves as an edit
recommender to fill knowledge gaps; and in web search, it offers lists of
salient statements about queried subjects for higher user engagement.

## Introduction

### Motivation

A community consists of a group of people who share a commonality such as
geography (Texans), religion (Christian), ethnicity (Arab), or a combination
(Arab Texans). One commonality that is less often discussed are communities of
passion or purpose, the so-called communities of interest (Fischer 2001). This
refers to groups of people who share a profession, practice, or an interest.
For instance, members of White House Coronavirus Task Force is a community of
practitioners in the medical field. Not to be confused with the much broader
community of all medical practitioners, we focus on contextualized groups of
people. In this case, people who were appointed by the White House for a
specific task. A second example is recipients of ACM Fellowship, rather than
all computer scientists. This allows for fine-grained analyses over various
topics, cultures, time frames, etc.

One standard task for understanding communities is identifying their
demographic factors. Demographics are statistical information about a
community that includes such factors as gender, occupation, linguistic
background, nationality, and location (Ashraf 2020). In geo-based communities
for example, identifying demographics can contribute in local policy making or
in understanding consumer behavior for national businesses. In communities of
interest, it could contribute to identifying under-represented groups or
understanding cultural differences between similar communities across
countries or continents. Moreover, compiling top demographics facilitates the
task of finding outliers, i.e., members that have different characteristics
than the majority, e.g., Deborah Birx is a female while 86% of the White House
Coronavirus Task Force members are male. These can contribute in studies of
under-represented groups in different settings.

While extensive research has been conducted on demographics of large geo-based
communities (Chambers 2020; Brass 1996) and topic-specific study cases
(Poschmann and Goldenstein 2022; Sun and Peng 2021; Zhou et al. 2020), to the
best of our knowledge, this is the first work to address communities of
interest, releasing data which covers 16 topics from 4 main domains, namely
culture, geography, history & society, and STEM.

Community | Demographics | Outliers
---|---|---
members of White House Coronavirus Task Force |  male republican |  Deborah Birx (female) Jerome Adams (independent)
members of Indian National Science Academy |  male indian |  Jörg Hacker (german)
winners of ACM Fellowship |  male american computer science |  Susan Nycum (female, lawyer) Calvin Gotlieb (canadian)

Table 1: Demographics and outliers of 3 communities.

### Contributions

We construct a dataset to capture demographics and outliers about communities
of interest from Wikidata (Vrandečić and Krötzsch 2014). Table 1 shows three
sample communities of interest, including winners of ACM Fellowship. While
most winners are male American computer scientists, two outliers are Susan
Nycum who is not male neither a computer scientist, but a female lawyer, and
Calvin Gotlieb who is not American but Canadian. In a nutshell: we collect
communities of interest from Wikidata using pre-defined properties
111https://www.wikidata.org/wiki/Wikidata:List_of_properties indicating
interest, inspired by the idea for collecting peer groups in (Arnaout et al.
2021). For example position held indicates people holding the same public
office and member of indicates people belonging to the same organization, club
or group. We then query its topic and domain using the Wiki-topic tool. In the
case of ACM Fellows, the domain is STEM and the topic is Computing. Next, we
query Wikidata for demographic statements about members of the community using
demography-describing properties in Wikidata, e.g., gender (P21). We compute
the most frequent value for each of the factors, e.g., occupation-computer
scientist (0.93) and gender-male (0.80). Finally, we use common factors to
identify outliers in the community. These are members whose characteristics do
not fully match the community’s demography profile. For example, members who
are not computer scientists. This results in a group-centric dataset. In
addition, we construct a subject-centric dataset, by merging for every
subject, the list of salient statements across multiple communities they are a
member of.

Our datasets can be downloaded from the Zenodo sharing service at:
https://doi.org/10.5281/zenodo.7410436. We release the data in an easy to
parse JSON format. We also publish a web interface for friendly browsing:
https://wikiknowledge.onrender.com/demographics.

On usage, we believe this dataset can provide useful seeds in social computing
problems, by identifying demographic data for a better analysis of communities
of interest and action-warranting under-representations, e.g., gender or
ethnicity-bias in certain communities. Downstream applications include
comparing similar communities within the same culture/country, e.g., the
majority of winners of the Latin Grammy Award for Best New Artist are men
while BET Award for Best New Artist are women (both American musical awards
for newcomers). Moreover, the data can give insights into equivalent
communities of interest between different cultures, e.g., the typical Prime
Minister of Iran is Shia but the typical Prime Ministers of Lebanon is Sunni,
which are the two main religious branches of Islam. One might be intrigued to
look more into the historical reasons of these findings or predict future
election winners. Beside social sciences, this data can be useful for
knowledge curation of collaborative knowledge bases such as Wikidata by
providing the curators with edit recommendations about missing information,
e.g., a missing profession of a member in a community can be derived from its
common profession (the typical Turing Award recipient is a computer
scientist). Finally, producing salient provenance-extended statements about
outliers can increase user engagement in web search (Tsurel et al. 2017),
e.g., unlike 560 out of 636 recorded winners of Presidential Medal of Freedom,
Stephen Hawking is not an American, but British, recipient.

## Dataset Creation

### Identifying Communities

We choose Wikidata as our source of communities of interest. It is a web-scale
collaborative knowledge base which covers 97% of Wikipedia articles (called
items). For example, ACM Fellowship as an article
(https://en.wikipedia.org/wiki/ACM_Fellow) and
(https://www.wikidata.org/wiki/Q18748039) as an item.

Wikidata, and structured stores of data in general, are a good source for
constructing communities of interest, due to their well-canonicalized entities
and properties, i.e., using IDs. For instance, one does not have to worry if
an Oscar is referred to using the word Oscar or Academy Award or different
wording but can simply query it using its unique ID Q19020. Moreover, Wikidata
contains additional information about every item, including temporal signals
and links to Wikipedia articles, which allows for more sophisticated use
cases, e.g., demographics in a certain era. We pick 6 properties indicating
interest or profession (position held \- P39, award \- P166, participant in \-
P1344, candidate in election \- P3602, nominated for \- P4353, member of \-
P463). We instantiate a SPARQL query 222https://query.wikidata.org/ with a
property of interest and one of its objects (community) to collect its
members, e.g., select distinct ?subject where {?subject wdt:P166 wd:Q18748039}
is used to collect ACM Fellowship recipients. A list of subjects is returned,
including Thomas Henzinger, Susan Nycum, Calvin Gotlieb, etc.

The outcome of this step is 7.5k communities of interest covering 16 topics
and 345k subjects. Given a community of interest, the Wiki-topic tool
333https://wiki-topic.toolforge.org/topic is queried for top-3 topics. An
overview is shown in Table 2 444For readability, we only display topics with
$\geq$ 50 communities.. Note that a community can belong to more than 1 topic,
e.g., Presidents of the Senate of Nigeria is both related to Politics under
History & Society and to Africa under Geography.

Topic | Domain | # of Communities | Sample Community (title, recorded members, sample member)
---|---|---|---
Biography |  Culture | 1358 |  Winners of Suffrage Science award, 35, Pippa Goldschmidt
Literature | 416 |  Winners of EU Prize for Literature, 96, Magdalena Parys
Media | 1071 |  Winners of Academy Award for Best Actress, 79, Emma Stone
Performing Arts | 74 |  Winners of Special Tony Award, 21, Judy Garland
Philosophy & Religion | 250 |  Popes, 269, John Paul II
Sports | 267 |  Winners of NBA Coach of the Year Award, 52, George Karl
Visual Arts | 189 |  Members of Royal Academy of Arts, 840, William Etty
Africa |  Geography | 191 |  Presidents of Uganda, 13, Idi Amin
Americas | 927 |  Governors of Wisconsin, 47, Jim Doyle
Asia | 947 |  Chiefs of the Philippine National Police, 23, Rodolfo Azurin Jr.
Europe | 2737 |  Prime Ministers of Poland, 44, Jan Olszewski
Oceania | 247 |  Winners of New Zealand Order of Merit, 148, Charles Higham
Business & Economics |  History | 148 |  Winners of Queen Elizabeth Prize for Engineering, 11, Tim Berners-Lee
Education | 134 |  Chancellors of the University of Oxford, 28, Roy Jenkins
History | 252 |  Dukes of Normandy, 26, William the Conqueror
Military & Warfare | 386 |  United States secretaries of War, 49, William Howard Taft
Politics & Government | 950 |  Leaders of the House of Commons, 80, Andrew Lansley
Society | 63 |  Winners of Civil Courage Prize, 36, Alexei Navalny
Biology |  STEM | 116 |  Winners of Darwin Medal, 67, Ernst Haeckel
Chemistry | 63 |  Winners of ACS Award in Pure Chemistry, 92, Stuart Schreiber
Earth & Environment | 76 |  Winners of Paleontological Society Medal, 54, Alfred Romer
Mathematics | 75 |  Winners of SIAM Fellow, 537, Ernst Hairer
Medicine & Health | 98 |  Winners of Pollin Prize for Pediatric Research, 17, Basil Hetzel
Physics | 120 |  Winners of Nobel Prize in Physics, 222, Albert Einstein
Space |  | 55 |  Winners of NASA Distinguished Service Medal, 319, Frank Borman

Table 2: Overview of covered topics and sample communities.

### Defining Demographic Factors

Now that we have communities of interest with their members and topics, we
want to identify their most frequent values given a set of standard
demographic factors (Ashraf 2020). We map each of those to equivalent Wikidata
properties (see Table 3). For instance, we identify the nationality of a
certain member using property P27.

### Inferring Demographics and Outliers

At this point, we have all the ingredients to start collecting community
demographics and outliers. For every community, e.g., award-ACM Fellow:

  1. 1.

From Wikidata, query values for the predefined demography-properties, e.g.,
gender(Thomas Henzinger, Calvin Gotlieb, ..) = male.

  2. 2.

Compute relative incidence of each factor-value pair, e.g., # male recipients
of ACM Fellowship/# recipients of ACM Fellowship= 673/839 = 0.80

  3. 3.

Sort by descending order of relative incidence, e.g., occupation-computer
scientist (0.93), gender-male (0.80), nationality-U.S. (0.58), etc..

  4. 4.

Collect outliers as members with demographic data not matching that of the
top-k of the community, e.g., NOT(gender-male) applies to Susan Nycum and
NOT(nationality-U.S.) applies to Calvin Gotlieb.

Accuracy of Inferred Information. When inferring non-asserted factors for
certain members, one unavoidable challenge is the correctness of these
inferences. Due to the open-world assumption knowledge bases such as Wikidata
postulates, absent statements can be either false (negative), or true but
simply absent (missing positive). Present statements in some cases can also be
undetectable using exact-match querying due to potential modelling issues. We
remedy these using three heuristics: (i) The partial completeness assumption
PCA (Galárraga et al. 2013; Dong et al. 2014), which asserts that if a subject
has at least one object for a given property, then there are no other objects
beyond those that are in the knowledge base, e.g., if we have at least 1 award
for subject X then we assume that their list of awards is complete. (ii)
Hierarchical checks, where we exploit the type system, i.e., class taxonomy,
in search for a contradiction of a certain negated factor. For instance,
occupation-Catholic priest does not hold for subject X, but occupation-Latin
Catholic priest does, and (Latin Catholic priest, subclassOf, Catholic priest)
is a statement in Wikidata. Hence, if X is a Latin Catholic priest, they are
also a Catholic priest. (iii) Semantic similarity checks to avoid possible
synonymous or near-synonymous contradictions, we compute the sentence
similarity between a candidate statement and an existing statement for subject
X. We do so using SBert (Reimers and Gurevych 2019) with 0.6 as a similarity
threshold, e.g., similarity (“teacher”, “professor”) = 0.62, avoiding the
inference that someone is a professor but not a teacher and vice versa. 95% of
the eliminated candidates are due to PCA, 2% due to hierarchical checks, and
due to 3% semantic similarity checks.

Factor | Wikidata Property
---|---
Gender | sex or gender (P21)
Sexual orientation | sexual orientation (P91)
Occupation | occupation (P106)
Political leaning | member of political party (P102)
Religion | religion or worldview (P140)
Linguistic background | native language (P103)
Ethnicity & race | ethnic group (P172)
Nationality | country of citizenship (P27)
Location | residence (P551)
Table 3: Standard demographic factors and their properties.

## Dataset Description

We release two datasets on Zenodo 555https://doi.org/10.5281/zenodo.7410436
and a browsing interface 666https://wikiknowledge.onrender.com/demographics/.

### Group-centric Dataset

This dataset consists of 7530 rows in English language, with a total size of
64MB in JSON format. The fields of a JSON record are:

  * •

Title ID: title of the community using Wikidata IDs.

  * •

Title label: title of the community using equivalent Wikidata labels.

  * •

Number of recorded members: number of subjects in Wikidata that belong to the
community.

  * •

Topics: a list of topics describing the community, e.g., Culture.Media.Music.

  * •

Demographic factors: a list of top demographics, each consisting of an ID,
value, and a score. The ID describes the factor using Wikidata identifiers and
value has their labels. The score is the relative incidence within the
community.

  * •

Outliers:

    * –

Reason: a statement on why the following members are considered outliers.

    * –

Score: a numerical value indicating the frequency of this factor in the
community.

    * –

Members: a list of members for which this factor does not hold.

A sample record 777For readability we omit some of the listed fields. from the
group-centric dataset:

⬇

1{

2 "title": "holders of position Lord Mayor of Dublin",

3 "recorded_members": 91,

4 "topics": ["Geography.Northern_Europe"],

5 "demographics": [

6 "gender-male",

7 "occupation-politician"

8 ],

9 "outliers": [

10 {

11 "reason": "NOT(male) unlike 81 out of 91 recorded members",

12 "members": [

13 "Catherine Byrne (female)",

14 "Emer Costello (female)",

15 "Alison Gilliland (female)"

16 ]

17 },

18 {

19 "reason": "NOT(politician) unlike 47 out of 91 recorded members",

20 "members": [

21 "John D’Arcy (businessperson)",

22 "Dermot Lacey (environmentalist)"

23 ]

24 }

25 ]

26}

Listing 1: JSON record from Group-centric dataset.

### Subject-centric Dataset

We rerun the method described in the previous section but with two
adjustments: for a given subject, we merge all outlier statements across
different communities and rank them by descending order of incidence.
Moreover, we extend the list of demography-properties to all possible Wikidata
properties. This subject-centric dataset consists of 345435 rows in English
language, with a total size of 172MB in JSON format.

The fields of a JSON record are:

  * •

Subject ID: the Wikidata ID of the subject.

  * •

Subject label: its equivalent label.

  * •

Statements: a list of salient statements across all communities of interest
this individual is a member of:

    * –

Statement ID: a statement using Wikidata ids.

    * –

Statement label: a statement using Wikidata labels.

    * –

Score: relative incidence.

A sample record 888For readability we omit some of the described fields. from
the subject-centric dataset:

⬇

1{

2 "subject": "Serena Williams",

3 "statements": [

4 {

5 "statement": "NOT(gender-male) but (female) 56 out of 68 recorded winners of
L’Equipe Champion of Champions.",

6 "score": 0.82

7 },

8 {

9 "statement": "NOT(sport-basketball) but (tennis) unlike 4 out of 8 recorded
winners of Best Female Athlete ESPY Award.",

10 "score": 0.50

11 }

12 ]

13}

Listing 2: JSON record from the Subject-centric dataset.

### Browsing Interface

For users who wish to sample rows from each topic, we publish a web interface
that can be accessed at: https://wikiknowledge.onrender.com/demographics/. A
screenshot is shown in Figure 1. It shows that the user selected Computing as
their topic of interest. It contains 40 communities about computing, one of
which is Winners of Turing Award with top demographic factors as male computer
scientists from the U.S.. Outliers include women, e.g., Barbara Liskov and
non-Americans, e.g., Tony Hoare who is from the U.K.. Users can also directly
search for communities of their choice. On top of this feature, the website
offers an entity summarization interface where users can query for the
favorite public figures and get a list of top salient statements about them.

Figure 1: Screenshot of the web interface, searching for demographics and
outliers of communities under Computing.

## Applications

### Demographic Data Analysis for Social Sciences

Discovering Under-represented Groups. One use case for our data is in academic
research in humanities. One standard social problem is identifying under-
represented groups (Atkinson et al. 2019; Burek 2021; Clark 1991). These
groups can be defined using one or more demographic factor, e.g., ethnicity,
gender. Our data can be considered a resource to answer questions such as: is
group X under-represented in community/domain Y? or what is the difference in
representation of group X between different domains? A more specific example
is shown in Figure 2. We show the fraction of female award recipients (in
Wikidata assigned gender-female) in STEM
999https://en.wikipedia.org/wiki/Wikipedia:AfC_sorting/STEM and in political
offices holders in different continents. These numbers can support needed
initiatives for more representation of certain groups, e.g., programs such as
Women in Tech. In this example, we used Wiki-topics to specify spatial and
topical dimensions, e.g., STEM.Physics for Physics awards and an intersection
of Regions.Europe and History_and_Society.Politics_and_government for
political offices in certain geographical regions.

Beyond the topics tool, our dataset is not isolated in terms of what we know
about communities of interest and their members. Each subject or group can be
linked to its Wikidata profile or Wikipedia article, allowing for more
customized analyses. For instance, in award winning or holding public offices,
statements are normally associated with temporal data in Wikidata, allowing
the user of this dataset to explore progress across time. For instance, the
charts in Figure 2 can be re-plotted to include time windows, i.e., female
award recipients in Physics [1960-1990], Physics [1991-present] and so on.

Figure 2: Female award recipients in STEM (left) and female political office
holders in different continents (right).

Exploring Cultural Differences in Governing. Our data can be used in political
science research such as understanding governing in different parts of the
world. One angle is to understand what kind of professions dominate public
offices, e.g., presidents, mayors, governors, ministers, … We compute these as
communities of interest created using the property position held
101010https://www.wikidata.org/wiki/Property:P39. We retain communities of
interest under Politics_and_Government, and assign each to its equivalent
geographical area, e.g., Central America (see Table 4). Note that we drop the
profession at rank 1, since it is politician for all the areas, and that
holding a certain office automatically turn one into a politician. This data
can give better understanding of certain cultures or in case of democracies,
how do people vote. Lawyer is a recurring profession in Americas, especially
in North America with quarter of public office holders with occupation-lawyer.

Area | Top professions
---|---
Central Africa | diplomat (0.27), economist (0.04), civil servant (0.01), philosopher (0.01), minister (0.01)
Eastern Africa | diplomat (0.09), judge (0.03), lawyer (0.03), military personnel (0.03), economist (0.02)
Northern Africa | diplomat (0.12), ruler (0.12), lawyer (0.03), military personnel (0.01), imam (0.01)
Southern Africa | judge (0.28), lawyer (0.11), civil servant (0.01), businessperson (0.01)
Western Africa | diplomat (0.17), lawyer (0.03), military personnel (0.03), economist (0.01), judge (0.01)
Central America | lawyer (0.07), diplomat (0.07), writer (0.02), economist (0.01), military personnel (0.01)
North America | lawyer (0.25), diplomat (0.06), judge (0.03), military personnel (0.01), businessperson (0.01)
South America | lawyer (0.17), diplomat (0.05), military personnel (0.02), journalist (0.01), historian (0.01)
East Asia | monarch (0.09), diplomat (0.07), lawyer (0.06), judge (0.06), prosecutor (0.01)
South Asia | diplomat (0.05), lawyer (0.03), economist (0.02), civil servant (0.02), judge (0.01)
Southeast Asia | sovereign (0.09), judge (0.08), lawyer (0.07), military personnel (0.03), diplomat (0.02)
West Asia | diplomat (0.12), sovereign (0.08), military personnel (0.05), physician (0.02), poet (0.01)
Eastern Europe | diplomat (0.12), economist (0.04), lawyer (0.02), monarch (0.02), university teacher (0.01)
Northern Europe | judge (0.08), diplomat (0.04), monarch (0.02), lawyer (0.02), journalist (0.01)
Southern Europe | diplomat (0.07), lawyer (0.04), military personnel (0.02), jurist (0.01), monarch (0.01)
Western Europe | lawyer (0.13), judge (0.06), diplomat (0.03), military personnel (0.02), suffragist (0.02), teacher (0.01)
Oceania | lawyer (0.08), diplomat (0.04), judge (0.01), pastoralist (0.01), solicitor (0.01), farmer (0.01)
Table 4: Top professions in political offices in different parts of the world.

### Edit Recommendations for Collaborative Encyclopedias

The Web-scale collaborative knowledge base Wikidata contains more than 100
million items (or subjects) which have received almost 2 billion edits since
its inception. Editors often need to prioritize their efforts, so useful tools
to guide them can improve data quality and completeness, e.g., the Recoin
plugin (Balaraman, Razniewski, and Nutt 2018) helps focus the editing on
missing properties of subjects. Our approach and dataset can be used to
improve this service by not only proposing relevant missing properties but
also proposing a full statement about that property. As mentioned, we consider
the PCA prior to inferring the negativity of a certain demographic factor. In
that step, one cannot assert absent information but can offer a calculated
guess of what that might be, leaving it for human curators to confirm or deny.
For example, Maja Vuković is a member of winners of IBM fellowship. For the
property occupation in Wikidata, she has zero values
111111https://www.wikidata.org/wiki/Q111536437 and Recoin lists occupation as
the top missing property. Given the demographic data we have about professions
of this community she is member of, we propose computer scientist,
mathematician, and engineer as top 3 candidates. This can especially
contribute to the completeness of information about long tail entities.
Moreover, for subjects who are members of multiple Communities of Interest,
one can merge similar demographic factors across communities and average
confidence scores.

### Entity Summarization

Web search results of queries about public figures can be improved by
including salient and sometimes surprising facts. It increases user engagement
(Tsurel et al. 2017) to augment question answering and entity summarization
results with did you know-like statements. This is where our second dataset,
i.e., the subject-centric, shines. A user can use the autocompletion field in
our web interface to query statements about subject of their choice. These are
often surprising or unexpected statements. Examples:

  * •

Did you know that unlike 96% of Oscar for Best Director winners, Bong Joon-ho
does not speak English, but Korean? 121212In fact he was accompanied by a
translator to deliver his acceptance speech.

  * •

Did you know that unlike 88% of winners of the Presidential Medal of Freedom,
Stephen Hawking is a non-American, but British, winner?

  * •

Did you know that unlike 91% of recipients of Liebig Medal (established by
Association of German Chemists), recipient Max Planck was not a chemist but a
physicist?

## Ethical & FAIR Considerations

In this work, datasets released are based on public information and tools
provided by Wikimedia projects. We do not use any personal information. Every
dataset record, i.e., JSON object, include all equivalent labels and IDs of
properties and items retrieved from Wikidata (December 2022 version). Users
can refer to these IDs for more details and definitions.

## Conclusion

In this work, we generate demographics from Wiki-based sources by exploiting
information about communities of interest. We release two datasets and publish
a web interface for friendly browsing. Finally, we show three purposes for the
data, namely as a resource for social sciences problems, edit recommender for
collaborative knowledge bases, and salient fact generator for a better search
experience.

## Acknowledgments

Funded by the German Research Foundation (DFG: Deutsche
Forschungsgemeinschaft) - Project 453095897.

## References

  * Arnaout et al. (2021) Arnaout, H.; Razniewski, S.; Weikum, G.; and Pan, J. Z. 2021.  Negative Statements Considered Useful.  _JWS_.
  * Ashraf (2020) Ashraf, M. A. 2020.  Demographic factors, compensation, job satisfaction and organizational commitment in private university: an analysis using SEM.  _Journal of Global Responsibility_.
  * Atkinson et al. (2019) Atkinson, R.; Lu, P.; Cho, N. L.; Melnitchouk, N.; and Kuo, L. E. 2019.  Gender disparities in award recipients from surgical specialty societies.  _Surgery_.
  * Balaraman, Razniewski, and Nutt (2018) Balaraman, V.; Razniewski, S.; and Nutt, W. 2018.  Recoin: Relative Completeness in Wikidata.  In _WWW Companion_. International World Wide Web Conferences Steering Committee.
  * Brass (1996) Brass, W. 1996.  Demographic Data Analysis in Less Developed Countries: 1946–1996.  _Population Studies_.
  * Burek (2021) Burek, C. V. 2021.  Female medal and fund recipients of the Geological Society of London: a historical perspective.  _Geological Society, London, Special Publications_.
  * Chambers (2020) Chambers, J. 2020.  Hybrid gridded demographic data for the world, 1950-2020.
  * Clark (1991) Clark, J. 1991.  Getting There: Women in Political Office.  _The ANNALS of the American Academy of Political and Social Science_.
  * Dong et al. (2014) Dong, X. L.; Gabrilovich, E.; Heitz, G.; Horn, W.; Murphy, K.; Sun, S.; and Zhang, W. 2014.  From Data Fusion to Knowledge Fusion.  _VLDB_.
  * Fischer (2001) Fischer, G. 2001.  Communities of interest: Learning through the interaction of multiple knowledge systems.  In _IRIS_.
  * Galárraga et al. (2013) Galárraga, L. A.; Teflioudi, C.; Hose, K.; and Suchanek, F. 2013.  AMIE: association rule mining under incomplete evidence in ontological knowledge bases.  In _WWW_.
  * Poschmann and Goldenstein (2022) Poschmann, P.; and Goldenstein, J. 2022.  Disambiguating and Specifying Social Actors in Big Data: Using Wikipedia as a Data Source for Demographic Information.  _Sociological Methods & Research_.
  * Reimers and Gurevych (2019) Reimers, N.; and Gurevych, I. 2019.  Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks.  In _EMNLP_.
  * Sun and Peng (2021) Sun, J.; and Peng, N. 2021.  Men are elected, women are married: Events gender bias on wikipedia.  _arXiv preprint arXiv:2106.01601_.
  * Tsurel et al. (2017) Tsurel, D.; Pelleg, D.; Guy, I.; and Shahaf, D. 2017.  Fun Facts: Automatic Trivia Fact Extraction from Wikipedia.
  * Vrandečić and Krötzsch (2014) Vrandečić, D.; and Krötzsch, M. 2014.  Wikidata: a Free Collaborative Knowledge base.  _CACM_.
  * Zhou et al. (2020) Zhou, S.-J.; Zhang, L.-G.; Wang, L.-L.; Guo, Z.-C.; Wang, J.-Q.; Chen, J.-C.; Liu, M.; Chen, X.; and Chen, J.-X. 2020.  Prevalence and socio-demographic correlates of psychological health problems in Chinese adolescents during the outbreak of COVID-19.  _European child & adolescent psychiatry_.


# Coherent-State Ladder Time-Dependent Variational Principle
for Open Quantum Systems

David S. Schlegel  0000-0003-2013-1676<EMAIL_ADDRESS>Laboratory of
Theoretical Physics of Nanosystems (LTPN), Institute of Physics, École
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Center for Quantum Science and Engineering, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Fabrizio Minganti
0000-0003-4850-1130 Laboratory of Theoretical Physics of Nanosystems (LTPN),
Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Center for Quantum Science and Engineering, École
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Vincenzo Savona  0000-0002-8984-6584 Laboratory of Theoretical Physics of
Nanosystems (LTPN), Institute of Physics, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Center for Quantum Science and
Engineering, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland

###### Abstract

We present a new paradigm for the dynamical simulation of interacting many-
boson open quantum systems. The method relies on a variational ansatz for the
$n$-boson density matrix, in terms of a superposition of photon-added coherent
states. It is most efficient for the simulation of driven-dissipative systems
where the state is well described by quantum fluctuations on top of a
displaced field, making it suitable for the simulation of several coupled
modes with large occupation numbers, that are otherwise very challenging using
a Fock-space expansion. We test our method on several examples, demonstrating
its potential application to the predictive simulation of interacting bosonic
systems and cat qubits.

#### Introduction.—

Quantum technologies rely on the generation and control of many-body quantum
states on platforms that maintain long-lived quantum coherence [1, 2, 3].
While quantum simulators and computers are expected to eventually outperform
classical emulation [4, 5, 6, 7], the certification and verification of
quantum protocols still require efficient approximate methods to simulate
entangled many-body dynamics [8, 9]. This task is especially challenging for
the many bosonic degrees of freedom characterizing superconducting circuits
[10, 11], integrated photonic platforms [12, 13], optomechanical resonators
[14], and trapped ions [15] among others.

In quantum optics, several methods have been developed to model weakly-
interacting systems close to their classical limit, based on the Gross-
Pitaevskii equation (GPE) [16, 17, 18] or the Gutzwiller mean-field
approximation [19, 20, 21, 22, 23]. Interactions, however, can induce quantum
correlations that cannot be captured by a mean-field approximation, nor by
methods such as the truncated Wigner approximation or low-order cumulant
expansion [24, 25]. Similarly, driven bosonic systems can operate in a regime
where the state inherits a large degree of coherence from the engineered
driving fields. A prototypical example are bosonic quantum codes, where drive
and dissipation create and stabilize quantum superpositions of coherent states
[26, 27, 28, 29]. The prohibitively large Hilbert space of multimode, driven
bosonic systems is often an obstacle to the use of conventional simulation
techniques based on Fock-space truncation [30, 31]. Novel methods, such as the
_shifted Fock basis_ [32], mitigate this problem when simulating the steady
state of the system subject to a stationary driving field. The efficient
numerical simulation of the dynamics of driven-dissipative many-body bosonic
systems remains, however, a major challenge.

In this Letter, we model the dynamics of a many-body driven-dissipative
bosonic system using the time-dependent variational principle [33, 34, 35]
(TDVP). We introduce an ansatz based on a set of unnormalized coherent states
called a coherent-state ladder. These states are chosen so to have optimal
overlap with the corner of the Hilbert space spanned by the system dynamics.

The method naturally accounts for symmetries, such as rotational symmetries in
phase space characterizing, e.g., bosonic cat codes [36, 26, 37, 38, 39, 32,
40] and other rotation symmetric bosonic codes [41, 42], thus permitting an
efficient study of time-dependent gate protocols and the occurrence of errors.
The number of basis states enters as a controllable parameter that sets the
accuracy of the approximation. It allows simulating coupled interacting
bosonic systems of considerable size requiring only moderate computational
resources.

Hereunder, we detail the principle of the method and benchmark its performance
by applying it to various models of driven-dissipative bosonic systems. Our
method outperforms the Fock-space truncation method in all examples, requiring
fewer basis states to represent the system’s dynamics accurately.

(a)

method |  large $\hat{n}$ |  inter-actions |  symme-tries |  dyna-mics
---|---|---|---|---
Fock-space [30, 31] |  no |  yes |  yes |  yes
GPE [16, 17] |  yes |  no |  no |  yes
shifted Fock [32] |  yes |  yes |  yes |  no
our method |  yes |  yes |  yes |  yes
(b)

Figure 1: (a) Schematic illustration of the variational ansatz. The TDVP
equations of motion in Eqs. (5-6) couple the coherent-field amplitude $\alpha$
(left) of the Bargmann state $\lVert{\alpha}\rangle$ [c.f. Eq. (3)] and the
density matrix elements $B_{ij}$ of the co-moving basis
$\mathcal{S}_{N}=\\{\lVert{\alpha}\rangle,\hat{a}^{\dagger}\lVert{\alpha}\rangle,{{}\hat{a}^{\dagger}}^{2}\lVert{\alpha}\rangle\dots{}\hat{a}^{\dagger}{}^{N}\lVert{\alpha}\rangle\\}$.
(b) Comparison of various methods with respect to the ability to capture large
occupation number $n$, interactions, rotational symmetries, and the dynamics
of an open bosonic quantum system.

#### Method.—

We consider a bosonic open quantum system coupled to a Markovian bath [43],
whose dynamics is governed by the Lindblad master equation

$\derivative{\hat{\rho}}{t}=-i\left[\hat{H},\hat{\rho}\right]+\sum\limits_{j}\gamma_{j}\mathcal{D}\left[\hat{L}_{j}\right](\hat{\rho}).$
(1)

Here, $\hat{H}$ is the Hamiltonian and $\mathcal{D}[\hat{L}_{j}](\hat{\rho})$
are the dissipators, with jump-operators $\hat{L}_{j}$ and dissipation rates
$\gamma_{j}$, defined by
$\mathcal{D}[\hat{A}](\hat{\rho})\equiv\hat{A}\hat{\rho}\hat{A}^{\dagger}-\frac{1}{2}\hat{A}^{\dagger}\hat{A}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{A}^{\dagger}\hat{A}$.

To implement the TDVP, we adopt the formalism proposed in Ref. [44]. The
ansatz for the density matrix is then

$\hat{\rho}=\sum_{i,j}^{N}B_{ij}(t)\outerproduct{\phi_{i}(t)}{\phi_{j}(t)},$
(2)

where the states $\ket{\phi_{i}(t)}$ are not in general orthogonal, and are
expressed in terms of variational parameters. The matrix elements $B_{ij}$
also enters as additional variational parameters. Once a criterion for
determining the states $\ket{\phi_{i}(t)}$ is introduced, the number $N$ will
determine the expressive power of the ansatz and the accuracy of the
variational approach.

The parametrized states $\ket{\phi_{i}(t)}$ that we introduce are based on the
_Bargmann_ states $\lVert{\alpha}\rangle$ [45, 46, 47], i.e., non-normalized
coherent states defined by

$\lVert{\alpha}\rangle=e^{\alpha\hat{a}^{\dagger}}\ket{0}=\sum\limits_{n=0}^{\infty}\frac{\alpha^{n}}{\sqrt{n!}}\ket{n}=e^{\frac{|\alpha|^{2}}{2}}\ket{\alpha},$
(3)

with $\ket{n}$ being Fock states. Starting from $\lVert{\alpha}\rangle$, the
subsequent application of the ladder operator $\hat{a}^{\dagger}$ yields the
(non-orthonormal) photon-added coherent-state ladder basis
$\mathcal{S}_{N}=\\{\lVert{\alpha}\rangle,\hat{a}^{\dagger}\lVert{\alpha}\rangle,{{}\hat{a}^{\dagger}}^{2}\lVert{\alpha}\rangle\dots{}\hat{a}^{\dagger}{}^{N}\lVert{\alpha}\rangle\\}$
[48]. For $\lVert{\alpha=0}\rangle$, $\mathcal{S}_{N}$ is just the Fock
(number) basis. For $\alpha\neq 0$, $\mathcal{S}_{N}$ describes excitations of
the state $\lVert{\alpha}\rangle$.

The states $\mathcal{S}_{N}$ bring several advantages in the simulation of the
dynamics. First, they express efficiently the state of a driven system, which
is expected to be given by a coherent state plus small quantum fluctuations.
Furthermore, the Bargmann states obey:

$\hat{a}\lVert{\alpha}\rangle=\alpha\lVert{\alpha}\rangle,\quad\partialderivative{\lVert{\alpha}\rangle}{\alpha}=\hat{a}^{\dagger}\lVert{\alpha}\rangle.$
(4)

From Eq. (4) it follows that the tangent space of $\mathcal{S}_{N}$ lies in
$\mathcal{S}_{N+1}$. In TDVP, the tangent space contains the gradient
function, which dictates the updates of the variational parameters, whose
evaluation is therefore greatly simplified. Notice also that Bargmann states
are holomorphic functions of the coherent-field amplitude $\alpha$, which
further simplifies the TDVP formulation using McLachlan’s principle [49, 50].

As depicted schematically in Fig. 1, we propose to use the McLachlan’s
variational principle [35, 49, 51] on a multi-mode bosonic system, by
extending the density matrix ansatz in Eq. (2) to $M$ bosonic modes. Each mode
is expressed by its coherent-field amplitude $\alpha_{k}$ [44] and by the
related coherent-state ladder
$\mathcal{S}_{N}^{(k)}=\\{\hat{a}_{k}^{n}\lVert{\alpha}\rangle,n=0\dots N\\}$,
where $k=1,\,\ldots,\,M$. The TDVP translates then into coupled differential
equations for the variational parameters $\boldsymbol{B}=\\{B_{i\,j}\\}$ and
$\alpha_{k}$ (See Supplemental Material for details)

$\displaystyle\boldsymbol{\dot{B}}$
$\displaystyle=\boldsymbol{S}^{-1}\boldsymbol{L}\boldsymbol{S}^{-1}-\boldsymbol{S}^{-1}\boldsymbol{\tau}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{\tau}^{\dagger}\boldsymbol{S}^{-1},$
(5) $\displaystyle\dot{\alpha_{k}}$
$\displaystyle=\frac{\text{Tr}\\{\boldsymbol{Y_{0}}^{(k)}\boldsymbol{B}\\}}{\text{Tr}\\{\boldsymbol{C_{0}}^{(k)}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}},$
(6)

where the matrices $\boldsymbol{S}$, $\boldsymbol{L}$, $\boldsymbol{C}$,
$\boldsymbol{\tau}$, $\boldsymbol{Y_{0}}^{(k)}$ and $\boldsymbol{C_{0}}^{(k)}$
depend on $\alpha_{k}$, $\dot{\alpha}_{k}$, and the model parameters.

#### Application Ia: Driven-dissipative Kerr resonator.—

As a first proof of principle, we consider a single system consisting of a
single-photon drive and Kerr-nonlinearity, whose Hamiltonian reads

$\hat{H}=\frac{U}{2}{{}\hat{a}^{\dagger}}^{2}\hat{a}^{2}+F(\hat{a}+\hat{a}^{\dagger}),$
(7)

with single-photon driving amplitude $F$ and Kerr-nonlinearity $U$. The system
is subject to single-photon loss $\kappa\mathcal{D}[\hat{a}]$. In the
classical limit of $U\rightarrow 0$ and $F\rightarrow\infty$, with $F\sqrt{U}$
constant, the required number of Fock-states needed to accurately describe the
systems’ dynamics increases when decreasing $U$, as the average photon number
of the system diverges as $(F/U)^{2/3}$ [52]. In Fig. 2(a), we show the number
of required basis states needed to reach a fidelity $\mathcal{F}>99\%$, both
for the Fock basis and our TDVP basis. Our results show that while the number
of Fock states increases by lowering $U$, the required number of basis states
for our TDVP method remains nearly constant. In Fig. 2(b), we show the
infidelity $1-\mathcal{F}$ as a function of time for different basis sizes
$N$, starting from the steady-state predicted by the semiclassical
approximation. The results illustrate how the quality of the TDVP solution
increases with $N$.

#### Application Ib: Asymmetrically driven nonlinear photonic dimer.—

Figure 2:  (a) Number of basis states $N$ needed to keep the fidelity above
the value $\mathcal{F}=0.99$ along the dynamics of the dissipative Kerr
resonator for the Fock-basis (blue) and the coherent-ladder method (orange)
for $F=1.5\sqrt{\kappa^{3}/U}$. (b-c) Infidelity $1-\mathcal{F}$ as a function
of dimensionless time $\kappa t$, for (b) different numbers of basis states
$N$ in the dissipative Kerr resonator and (c) in the asymmetrically driven
nonlinear photonic dimer for $J=1.2\kappa$, $\Delta_{1,2}=2\kappa$,
$\alpha_{1}(t=0)=-1.0-1.84i$, $\alpha_{2}(t=0)=1.36+0.8i$.

Second, we consider a system of two coupled modes with Hamiltonian

$\displaystyle\hat{H}=$
$\displaystyle\sum_{i=1,2}-\Delta\hat{a}_{i}^{\dagger}\hat{a}_{i}+\frac{U}{2}{{}\hat{a}_{i}^{\dagger}}^{2}\hat{a}_{i}^{2}$
(8)
$\displaystyle-J\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{1}\hat{a}_{2}^{\dagger}\right)+F\left(\hat{a}_{1}^{\dagger}+\hat{a}_{1}\right),$

and single-photon loss processes $\kappa\mathcal{D}[\hat{a}_{1}]$ and
$\kappa\mathcal{D}[\hat{a}_{2}]$. Here, $\Delta$ is the detuning of the
resonators with respect to the drive frequency, $J$ the hopping interaction,
and $F$ the driving-field amplitude, acting only on one of the two modes. A
semiclassical treatment using the GPE predicts multiple parametrically
unstable regions for different parameter regimes [53, 54]. Again, in the
classical limit of $U\to 0$ and $F\to\infty$, with $F\sqrt{U}$ constant,
accurately describing the effect of quantum fluctuations with a truncated-Fock
approach becomes challenging as $|\alpha|^{2}$, and thus the average number of
bosons per mode, diverge as $(F/U)^{2/3}$ [52]. Furthermore, the field
amplitude $\alpha$ varies periodically along the dynamics, and cannot be
modeled using a stationary ansatz, like e.g. with the shifted Fock basis [32].
This system, therefore, provides an ideal testbed for our method.

Starting from an initial coherent state, we simulate the time-evolution for
different basis sizes $N$. The results are shown in Fig. 2(c), where the
infidelity $1-\mathcal{F}$, with respect to the exact solution, is plotted.
Our results indicate that with few basis states in the variational ansatz, the
essential dynamics of the system can be accurately captured.

#### Rotational symmetries.—

Our approach conveniently facilitates the inclusion of rotational symmetries,
thereby enabling the modeling of rotation-symmetric quantum codes such as cat
codes. As an example, we consider a $\mathds{Z}_{2}$ rotational symmetry whose
symmetry subspaces can be decomposed in even and odd superpositions of
Bargmann states with opposite displacement, yielding the unnormalized even and
odd cat states

$\lVert{\mathcal{C}_{\alpha}^{\pm}}\rangle=\lVert{\alpha}\rangle\pm\lVert{-\alpha}\rangle.$
(9)

Similar to the previous case, we can construct a cat-ladder basis comprising
the two symmetry sectors by applying $\hat{a}^{\dagger}$ on the even and odd
cat states in Eq. (9). Note that
$\partialderivative{\lVert{\mathcal{C}_{\alpha}^{\pm}}\rangle}{\alpha}=\hat{a}^{\dagger}\lVert{\mathcal{C}_{\alpha}^{\mp}}\rangle$.
Equations (5) and (6) maintain the same structure, with slight modifications
in the involved matrices (See Supplemental Material for details). A
$\mathds{Z}_{N}$ rotational symmetry can similarly be encoded by defining a
basis of superpositions of rotated coherent Bargman states
$e^{ik2\pi\hat{n}/N}\lVert{\alpha}\rangle$, with $k=0,\,\dots,\,N-1$. In the
following, we apply our method to systems of interacting cat qubits.

#### Application IIa: Dynamics of a quenched cat qubit.—

Let us consider the dynamics of a cat qubit, stabilized by the simultaneous
action of two-photon dissipation and Kerr nonlinearity when the drive is
quenched. The system is described by the Hamiltonian

$\hat{H}=\left(\frac{G(t)}{2}\hat{a}^{2}+\frac{G(t)^{*}}{2}\hat{a}^{\dagger
2}\right)+\frac{U}{2}{{}\hat{a}^{\dagger}}^{2}\hat{a}^{2},$ (10)

and dissipator $\eta\mathcal{D}[\hat{a}^{2}]$, where $G(t)=G_{0}$ for $t<0$
and $G(t)=G_{1}$ for $t>0$ is the two-photon drive amplitude and $\eta$ the
two-photon dissipation rate. We suppose that the system reached its steady
state at $t<0$, and we investigate the dynamics after the quench at $t=0$.

Figure 3: (a) Snapshots of the Wigner function at different times along the
dynamics for an initial cat state, highlighting leakage out of the cat-state
manifold during the transient dynamics for $\eta=U/4$, $G_{0}=5iU$,
$G_{1}=-i5U$. (b) Infidelity $1-\mathcal{F}$ as a function of time for
different basis sizes $N$. (c) Parity expectation value
$\expectationvalue*{\hat{\Pi}}$ as a function of $t\Omega_{X}/\pi$ for
different ratios $F/\eta$, assuming an initial even cat state (solid lines)
and odd cat state (dashed lines), with $G=10\eta$. (d) For a fixed ratio
$F/\eta=4$, the quantity $\expectationvalue*{\hat{\Pi}}$ is plotted for
different basis sizes $N$.

In Fig. 3(a), we depict snapshots of the Wigner functions along the time
evolution. The leakage in the transient dynamics appears as a departure from
the shape of an ideal cat state. We show the effect of the basis cutoff $N$ in
Fig. 3(b). Despite the large leakage out of the cat-qubit manifold, a
relatively small basis size captures the whole dynamics with high accuracy.

#### Application IIb: Coherent non-adiabatic rotation of a cat qubit.—

We consider now the case $U=0$, where the steady-state is spanned by the cat-
qubit manifold with amplitude $\alpha=\sqrt{iG/{\eta}}$. We encode
$\ket{{0/1}_{\mathrm{L}}}\equiv\ket*{\mathcal{C}_{\alpha}^{\pm}}$ as the
logical $\hat{Z}$-eigenstates. These states are eigenstates of the photon
number parity $\hat{\Pi}=(-1)^{\hat{n}}$. A coherent rotation around the
logical $\hat{X}$-axis can be realized by a single-photon drive
$\hat{H}_{X}=F(\hat{a}+\hat{a}^{\dagger})$ [36], with single-photon driving
amplitude $F$. In the limit $F\ll\eta$, the code space remains unperturbed,
and the qubit is adiabatically rotated by an angle $\varphi=\Omega_{X}t$, with
Rabi frequency $\Omega_{X}=\sqrt{2}F|\alpha|$ [36]. To analyze the non-
adiabatic effect of the rotation-gate $\hat{H}_{X}$, we simulate the dynamics
of the system starting from $\ket*{\mathcal{C}_{\alpha}^{\pm}}$. In Fig. 3(c),
we show the average parity $\expectationvalue*{\hat{\Pi}}$ for different
single-photon driving amplitudes $F$ as a function of time. For $F/\eta\gtrsim
1$, the single-photon drive induces non-adiabatic processes, resulting in code
leakage and a reduction of the fidelity of the rotation gate. To assess the
efficiency of the method, we plot in Fig. 3(d) the parity
$\expectationvalue*{\hat{\Pi}}$ for different basis sizes for a fixed driving
strength $F$. Again, our results reveal that even for a moderate truncation,
the non-adiabatic nature of the $\hat{X}$-rotation gate can be adequately
described, significantly reducing the computational overhead compared to a
truncated Fock approach.

#### Application III: Interacting cat qubit dynamics.—

We now consider a two-mode cat-like bosonic system. Each mode ($\hat{a}_{1}$
and $\hat{a}_{2}$) is described by a two-photon driven-dissipative Kerr
resonator analogous to Eq. (10) that interacts through the hopping (beam-
splitter) Hamiltonian

$\hat{H}_{\mathrm{int}}=\sum\limits_{\langle
i,j\rangle}J_{ij}(\hat{a}_{i}\hat{a}_{j}^{\dagger}+\hat{a}_{i}^{\dagger}\hat{a}_{j}),$
(11)

with hopping interaction $J_{ij}$. The system conserves the total photon
number parity $\expectationvalue*{\hat{\Pi}_{1}\hat{\Pi}_{2}}$ [55]. The
coherent tunneling interaction $\hat{H}_{\mathrm{int}}$ enables coherent
exchange of excitations between the two modes, thereby encoding an entangling
two-qubit gate. We simulate the time evolution of the system starting from an
initial product state
$\ket{\psi(0)}=\ket*{\mathcal{C}_{\alpha_{1}}^{+}}\otimes\ket*{\mathcal{C}_{\alpha_{2}}^{+}}$,
where $\alpha_{1}$ and $\alpha_{2}$ differ from those characterizing the
steady state of the cat code. This may represent a system where the drive
changes upon application of the $\hat{X}_{1}\hat{X}_{2}$ gate. We compute the
time evolution of the system, assuming different basis sizes $N$, and compare
to the exact result in Fig. 4(a), where we show the evolution of
$\expectationvalue*{\hat{a}^{2}}$ and the single-mode parity
$\expectationvalue*{\hat{\Pi}_{1}}$ in Fig. 4(b). Convergence to the exact
result arises already for small basis size $N$ per mode.

Figure 4: (a) Real and imaginary parts of $\expectationvalue*{\hat{a}^{2}}$,
as computed along the dynamics of a system of two interacting driven-
dissipative Kerr-resonators. (b) Single-mode parity
$\expectationvalue*{\hat{\Pi}_{1}}$ as a function of time for the same system,
setting $G_{i}=5U,\;\eta_{i}=0.25U,\;J_{12}=U$ and initial state:
$\ket{\psi(0)}=\ket*{\mathcal{C}_{\alpha_{1}}^{+},\mathcal{C}_{\alpha_{2}}^{+}}$,
with $\alpha_{1}=\alpha_{2}=2$. (c)-(d) Single-mode parity
$\expectationvalue*{\hat{\Pi}_{1}}$ and two-mode parity
$\expectationvalue*{\hat{\Pi}_{1}\hat{\Pi}_{2}}$ as a function of time, as
computed for a system of three coupled modes, setting
$G_{i}=5U\;\eta_{i}=0.25U,\;J_{12}=J_{23}=0.8U$ and initial state
$\ket{\psi(0)}=\ket*{\mathcal{C}_{\alpha_{1}}^{+},\mathcal{C}_{\alpha_{2}}^{+},\mathcal{C}_{\alpha_{3}}^{+}}$,
with $\alpha_{i}=2$.

Finally, we simulate a system consisting of three interacting bosonic modes
arranged in an open chain configuration, where the nearest neighbor hopping
constants $J_{12}$ and $J_{23}$ are non-zero. In the regime we consider,
simulating this system with a truncated Fock approach starts being
computationally challenging. Convergence vs $N$ is then assessed by simply
comparing results for increasing $N$. Figs. 4(c-d) display the time evolution
of the single-mode parity $\expectationvalue*{\hat{\Pi}_{1}}$ and two-mode
parity $\expectationvalue*{\hat{\Pi}_{1}\hat{\Pi}_{2}}$ for various single-
mode basis sizes $N$. Analogously to the two-mode case, the total photon
number parity $\expectationvalue*{\hat{\Pi}_{1}\hat{\Pi}_{2}\hat{\Pi}_{3}}$ is
preserved throughout the system dynamics. Even for the three-mode case,
convergence is reached within a reasonable number of basis states,
highlighting the potential of the method.

#### Conclusions.—

We have introduced an approach for simulating the dynamics of interacting
many-boson driven-dissipative systems. The coherent-state ladder approach
employs a variational ansatz for the $n$-boson density matrix, where both the
matrix elements and the displacement fields serve as time-dependent
variational parameters. This approach bridges the gap between quantum and
classical realms, making it suitable for simulating systems with large
occupation numbers that are beyond the capabilities of Fock-space methods.

We have demonstrated the potential of the method through several examples,
showcasing its applicability to interacting bosonic systems and cat qubits. It
outperforms Fock-space methods in terms of the number of basis states needed
to accurately represent the system’s dynamics.

Our method can straightforwardly be extended to include different
parametrizations of the basis, such as squeezed states or higher-order
processes, as well as other types of symmetries (See Supplemental Material for
details). It can thus accurately study the system dynamics of multi-mode
bosonic systems in the presence of quantum correlations.

The coherent-state ladder thus offers new possibilities for accurately
studying system dynamics of a wide range of single- and multi-mode bosonic
systems, whose simulation would be otherwise very challenging. It opens new
pathways in the simulation of systems of many coupled bosonic qubits to
efficiently optimize gate and error-correction protocols [38, 39, 32, 40].

###### Acknowledgements.

#### Acknowledgments.—

This work was supported by the Swiss National Science Foundation through
Projects No. 200020_185015 and 200020_215172 and was conducted with the
financial support of the EPFL Science Seed Fund 2021.

## References

  * Schoelkopf and Girvin [2008] R. J. Schoelkopf and S. M. Girvin,  _Wiring up quantum systems_ , Nature 451, 664 (2008).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard and S. Nascimbène,  _Quantum simulations with ultracold quantum gases_ , Nature Physics 8, 267 (2012).
  * Monroe and Kim [2013] C. Monroe and J. Kim,  _Scaling the Ion Trap Quantum Processor_ , Science 339, 1164 (2013).
  * Cirac and Zoller [2012] J. I. Cirac and P. Zoller,  _Goals and opportunities in quantum simulation_ , Nature physics 8, 264 (2012).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner et al.,  _Probing many-body dynamics on a 51-atom quantum simulator_ , Nature 551, 579 (2017).
  * Tacchino  _et al._ [2020] F. Tacchino, A. Chiesa, S. Carretta and D. Gerace,  _Quantum Computers as Universal Quantum Simulators: State-of-the-Art and Perspectives_ , Advanced Quantum Technologies 3, 1900052 (2020).
  * Cirac [2021] J. I. Cirac,  _Quantum computing and simulation_ , Nanophotonics 10, 453 (2021).
  * Montanaro [2016] A. Montanaro,  _Quantum algorithms: an overview_ , npj Quantum Information 2, 1 (2016).
  * Preskill [2018] J. Preskill,  _Quantum Computing in the NISQ era and beyond_ , Quantum 2, 79 (2018).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell et al.,  _Quantum supremacy using a programmable superconducting processor_ , Nature 574, 505 (2019).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson and W. D. Oliver,  _Superconducting Qubits: Current State of Play_ , Annu. Rev. Condens. Matter Phys. 11, 369 (2020).
  * Wang  _et al._ [2020] J. Wang, F. Sciarrino, A. Laing and M. G. Thompson,  _Integrated photonic quantum technologies_ , Nature Photonics 14, 273 (2020).
  * Pelucchi  _et al._ [2022] E. Pelucchi, G. Fagas, I. Aharonovich, D. Englund, E. Figueroa, Q. Gong, H. Hannes, J. Liu, C.-Y. Lu, N. Matsuda, J.-W. Pan, F. Schreck, F. Sciarrino, C. Silberhorn, J. Wang and K. D. Jöns,  _The potential and global outlook of integrated photonics for quantum technologies_ , Nature Reviews Physics 4, 194 (2022).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg and F. Marquardt,  _Cavity optomechanics_ , Rev. Mod. Phys. 86, 1391 (2014).
  * Bruzewicz  _et al._ [2019] C. D. Bruzewicz, J. Chiaverini, R. McConnell and J. M. Sage,  _Trapped-ion quantum computing: Progress and challenges_ , Applied Physics Reviews 6, 021314 (2019).
  * Gross [1961] E. P. Gross,  _Structure of a quantized vortex in boson systems_ , Il Nuovo Cimento (1955-1965) 20, 454 (1961).
  * Pitaevskii [1961] L. P. Pitaevskii,  _Vortex lines in an imperfect Bose gas_ , Sov. Phys. JETP. 13, 451 (1961).
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti,  _Quantum fluids of light_ , Rev. Mod. Phys. 85, 299 (2013).
  * Le Boité  _et al._ [2013] A. Le Boité, G. Orso and C. Ciuti,  _Steady-State Phases and Tunneling-Induced Instabilities in the Driven Dissipative Bose-Hubbard Model_ , Phys. Rev. Lett. 110, 233601 (2013).
  * Lee  _et al._ [2013] T. E. Lee, S. Gopalakrishnan and M. D. Lukin,  _Unconventional Magnetism via Optical Pumping of Interacting Spin Systems_ , Phys. Rev. Lett. 110, 257204 (2013).
  * Jin  _et al._ [2016] J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio and D. Rossini,  _Cluster Mean-Field Approach to the Steady-State Phase Diagram of Dissipative Spin Systems_ , Phys. Rev. X 6, 031011 (2016).
  * Casteels  _et al._ [2018] W. Casteels, R. M. Wilson and M. Wouters,  _Gutzwiller Monte Carlo approach for a critical dissipative spin model_ , Phys. Rev. A 97, 062107 (2018).
  * Huybrechts and Wouters [2020] D. Huybrechts and M. Wouters,  _Dynamical hysteresis properties of the driven-dissipative Bose-Hubbard model with a Gutzwiller Monte Carlo approach_ , Phys. Rev. A 102, 053706 (2020).
  * Sánchez-Barquilla  _et al._ [2020] M. Sánchez-Barquilla, R. E. F. Silva and J. Feist,  _Cumulant expansion for the treatment of light–matter interactions in arbitrary material structures_ , J. Chem. Phys. 152, 034108 (2020).
  * Vicentini  _et al._ [2018] F. Vicentini, F. Minganti, R. Rota, G. Orso and C. Ciuti,  _Critical slowing down in driven-dissipative Bose-Hubbard lattices_ , Phys. Rev. A 97, 013853 (2018).
  * Ofek  _et al._ [2016] N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret and R. J. Schoelkopf,  _Extending the lifetime of a quantum bit with error correction in superconducting circuits_ , Nature (London) 536, 441 (2016).
  * Rosenblum  _et al._ [2018] S. Rosenblum, P. Reinhold, M. Mirrahimi, L. Jiang, L. Frunzio and R. J. Schoelkopf,  _Fault-tolerant detection of a quantum error_ , Science 361, 266 (2018).
  * Sivak  _et al._ [2023] V. V. Sivak, A. Eickbusch, B. Royer, S. Singh, I. Tsioutsios, S. Ganjam, A. Miano, B. L. Brock, A. Z. Ding, L. Frunzio, S. M. Girvin, R. J. Schoelkopf and M. H. Devoret,  _Real-time quantum error correction beyond break-even_ , Nature 616, 50 (2023).
  * Joshi  _et al._ [2021] A. Joshi, K. Noh and Y. Y. Gao,  _Quantum information processing with bosonic qubits in circuit QED_ , Quantum Science and Technology 6, 033001 (2021).
  * Fock [1932] V. Fock,  _Konfigurationsraum und zweite Quantelung_ , Zeitschrift für Physik 75, 622 (1932).
  * Reed and Simon [1975] M. Reed and B. Simon,  _Fourier Analysis, Self-Adjointness_, 1st ed., Vol. 2 (Elsevier Science & Technology, 1975).
  * Chamberland  _et al._ [2022] C. Chamberland, K. Noh, P. Arrangoiz-Arriola, E. T. Campbell, C. T. Hann, J. Iverson, H. Putterman, T. C. Bohdanowicz, S. T. Flammia, A. Keller, G. Refael, J. Preskill, L. Jiang, A. H. Safavi-Naeini, O. Painter and F. G. Brandão,  _Building a Fault-Tolerant Quantum Computer Using Concatenated Cat Codes_ , PRX Quantum 3, 010329 (2022).
  * Dirac [1982] P. A. M. Dirac,  _The Principles of Quantum Mechanics_ , 4th ed. (Clarendon Press, Oxford, 1982).
  * Frenkel [1934] J. Frenkel,  _Wave Mechanics: Advanced General Theory_ , oxford at the clarendon press; 1st edition ed. (Oxford University Press, 1934).
  * McLachlan [1964] A. McLachlan,  _A variational solution of the time-dependent Schrodinger equation_ , Molecular Physics 8, 39 (1964).
  * Mirrahimi  _et al._ [2014] M. Mirrahimi, M. Leghtas, V. Albert, S. Touzard, R. Schoelkopf, L. Jiang and M. Devoret,  _Dynamically protected cat-qubits: a new paradigm for universal quantum computation_ , New Journal of Physics 16, 045014 (2014).
  * Albert  _et al._ [2018] V. V. Albert, K. Noh, K. Duivenvoorden, D. J. Young, R. T. Brierley, P. Reinhold, C. Vuillot, L. Li, C. Shen, S. M. Girvin, B. M. Terhal and L. Jiang,  _Performance and structure of single-mode bosonic codes_ , Physical Review A 97, 032346 (2018).
  * Guillaud and Mirrahimi [2019] J. Guillaud and M. Mirrahimi,  _Repetition Cat Qubits for Fault-Tolerant Quantum Computation_ , Phys. Rev. X 9, 041053 (2019).
  * Puri  _et al._ [2020] S. Puri, L. St-Jean, J. A. Gross, A. Grimm, N. E. Frattini, P. S. Iyer, A. Krishna, S. Touzard, L. Jiang, A. Blais, S. T. Flammia and S. M. Girvin,  _Bias-preserving gates with stabilized cat qubits_ , Science Advances 6, eaay5901 (2020).
  * Gautier  _et al._ [2023] R. Gautier, M. Mirrahimi and A. Sarlette,  _Designing High-Fidelity Gates for Dissipative Cat Qubits_ (2023), arXiv:2303.00760 [quant-ph] .
  * Schlegel  _et al._ [2022] D. S. Schlegel, F. Minganti and V. Savona,  _Quantum error correction using squeezed Schrödinger cat states_ , Phys. Rev. A 106, 022431 (2022).
  * Grimsmo  _et al._ [2020] A. L. Grimsmo, J. Combes and B. Q. Baragiola,  _Quantum Computing with Rotation-Symmetric Bosonic Codes_ , Phys. Rev. X 10, 011058 (2020).
  * Breuer and Petruccione [2007] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Joubert-Doriol and Izmaylov [2015] L. Joubert-Doriol and A. F. Izmaylov,  _Problem-free time-dependent variational principle for open quantum systems_ , The Journal of Chemical Physics 142, 134107 (2015).
  * Bargmann [1961] V. Bargmann,  _On a Hilbert space of analytic functions and an associated integral transform part I_ , Communications on Pure and Applied Mathematics 14, 187 (1961).
  * Segal [1967] I. Segal,  _Mathematical Problems of Relativistic Physics: With an Appendix on Group Representations in Hilbert Space_ , Lectures in applied mathematics; proceedings of the Summer seminar, Boulder, Colorado, 1960, 2 (American Mathematical Society, 1967).
  * Pérez-Esteva and Villegas-Blas [2000] S. Pérez-Esteva and C. Villegas-Blas,  _First Summer School in Analysis and Mathematical Physics: quantization, the Segal-Bargmann transform, and semiclassical analysis_ (American Mathematical Society, 2000).
  * Agarwal and Tara [1991] G. S. Agarwal and K. Tara,  _Nonclassical properties of states generated by the excitations on a coherent state_ , Phys. Rev. A 43, 492 (1991).
  * Broeckhove  _et al._ [1988] J. Broeckhove, L. Lathouwers, E. Kesteloot and P. Van Leuven,  _On the equivalence of time-dependent variational principles_ , Chemical Physics Letters 149, 547 (1988).
  * Kan [1981] K.-K. Kan,  _Equivalence of time-dependent variational descriptions of quantum systems and Hamilton’s mechanics_ , Phys. Rev. A 24, 2831 (1981).
  * Raab and Meyer [2000] A. Raab and H.-D. Meyer,  _Multiconfigurational expansions of density operators: equations of motion and their properties_ , Theoretical Chemistry Accounts 104 (2000).
  * Bartolo  _et al._ [2016] N. Bartolo, F. Minganti, W. Casteels and C. Ciuti,  _Exact steady state of a Kerr resonator with one- and two-photon driving and dissipation: Controllable Wigner-function multimodality and dissipative phase transitions_ , Phys. Rev. A 94, 033841 (2016).
  * Sarchi  _et al._ [2008] D. Sarchi, I. Carusotto, M. Wouters and V. Savona,  _Coherent dynamics and parametric instabilities of microcavity polaritons in double-well systems_ , Phys. Rev. B 77, 125324 (2008).
  * Seibold  _et al._ [2020] K. Seibold, R. Rota and V. Savona,  _Dissipative time crystal in an asymmetric nonlinear photonic dimer_ , Phys. Rev. A 101, 033839 (2020).
  * Albert and Jiang [2014] V. V. Albert and L. Jiang,  _Symmetries and conserved quantities in Lindblad master equations_ , Phys. Rev. A 89, 022118 (2014).
  * Pantaleoni  _et al._ [2020] G. Pantaleoni, B. Q. Baragiola and N. C. Menicucci,  _Modular Bosonic Subsystem Codes_ , Phys. Rev. Lett. 125, 040501 (2020).

Supplemental Material:
Coherent-State Ladder Time-Dependent Variational Principle
for Open Quantum Systems

##  I Derivation of the TDVP equations of motion

In the McLachlan TDVP for density matrices, the equations of motion arise from
the minimization of the Frobenius norm of the error [35, 49]

$\left\|\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}]\right\|=\text{Tr}\\{(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])^{\dagger}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\\}^{\frac{1}{2}},$
(12)

which is equivalent to the stationary condition [51]

$\text{Tr}\\{\delta\hat{\rho}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\\}=0,$
(13)

with an infinitesimal variation in the density matrix $\delta\hat{\rho}$. For
our particular choice of ansatz in Eq. (2), infinitesimal changes in the
density matrix $\hat{\rho}$ translate to the variational paramters as

$\delta\hat{\rho}=\sum\limits_{ij}\partialderivative{\hat{\rho}}{B_{ij}}\delta
B_{ij}+\sum\limits_{k}\partialderivative{\hat{\rho}}{\alpha_{k}}\delta\alpha_{k}+\sum\limits_{k}\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}\delta\alpha_{k}^{*}.$
(14)

Inserting the above expression in Eq. (13), and due to the independent nature
of the derivatives, we have three stationary conditions:

$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{B_{ij}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
$\displaystyle=0$ (15)
$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
$\displaystyle=0$ (16)
$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
$\displaystyle=0$ (17)

Since the last two are complex conjugates of one another, they constitute the
same condition. We define the overlap matrix

$[\boldsymbol{S}]_{ij}=\innerproduct{\phi_{i}}{\phi_{j}}.$ (18)

As the coherent-state- and cat-ladder basis states are not orthogonal,
$\boldsymbol{S}$ is not an orthogonal matrix. Note that due to this non-
orthonormal basis, the projector $\hat{P}$, projecting onto the ansatz basis
$\mathcal{S}$, is given by

$\hat{P}=\sum\limits_{ij}\ket{\phi_{i}}[\boldsymbol{S}^{-1}]_{ij}\bra{\phi_{j}}.$
(19)

We can rewrite Eq. (14) as

$\derivative{\hat{\rho}}{t}=\sum\limits_{cd}\dot{B}_{cd}\outerproduct*{\phi_{c}}{\phi_{d}}+B_{cd}\outerproduct*{\dot{\phi_{c}}}{\phi_{d}}+B_{cd}\outerproduct*{\phi_{c}}{\dot{\phi_{d}}}$
(20)

Replacing this expression into Eq. (15) yields:

$\displaystyle\text{Tr}\left\\{\partialderivative{\hat{\rho}}{B_{ij}}\left(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}]\right)\right\\}=0$
(21)
$\displaystyle=\sum\limits_{kl}\bra{\phi_{l}}\boldsymbol{S}_{kl}^{-1}\left[\outerproduct{\phi_{i}}{\phi_{j}}\left(\sum\limits_{cd}\dot{B}_{cd}\outerproduct*{\phi_{c}}{\phi_{d}}\right.\right.$
(22)
$\displaystyle\qquad\left.\left.+B_{cd}\outerproduct*{\dot{\phi_{c}}}{\phi_{d}}+B_{cd}\outerproduct*{\phi_{c}}{\dot{\phi_{d}}}-\mathcal{L}[\hat{\rho}]\right)\right]\ket{\phi_{k}}$
$\displaystyle=\sum\limits_{cd}S_{jc}\dot{B}_{cd}S_{di}+\tau_{jc}B_{cd}S_{di}+S_{jc}B_{cd}\tau_{id}^{*}-L_{ji},$
(23)

where we have introduced

$\displaystyle[\boldsymbol{\tau}]_{mn}$
$\displaystyle=\sum\limits_{k}\innerproduct{\phi_{m}}{\partialderivative{\phi_{n}}{\alpha_{k}}}\dot{\alpha_{k}},$
(24) $\displaystyle[\boldsymbol{L}]_{ij}$
$\displaystyle=\matrixelement{\phi_{i}}{\mathcal{L}[\hat{\rho}]}{\phi_{j}}.$
(25)

We thus have the condition

$\boldsymbol{S}\boldsymbol{\dot{B}}\boldsymbol{S}+\boldsymbol{\tau}\boldsymbol{B}\boldsymbol{S}+\boldsymbol{S}\boldsymbol{B}\boldsymbol{\tau}^{\dagger}-\boldsymbol{L}=0,$
(26)

which gives the equation of motion for the density matrix coefficients,

$\boldsymbol{\dot{B}}=\boldsymbol{S}^{-1}\boldsymbol{L}\boldsymbol{S}^{-1}-\boldsymbol{S}^{-1}\boldsymbol{\tau}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{\tau}^{\dagger}\boldsymbol{S}^{-1}.$
(27)

For the equation of motion for $\dot{\alpha}_{k}$, we have by inserting the
definition of $\dot{\hat{\rho}}$ from Eq. (20),

$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}=0$
(28)
$\displaystyle\begin{split}&=\sum\limits_{ab}\bra{\partial_{\alpha_{k}}\phi_{b}}\left\\{\sum\limits_{cd}\left(\dot{B}_{cd}\outerproduct*{\phi_{c}}{\phi_{d}}\right.\right.\\\
&\qquad+\left.\left.B_{cd}\outerproduct*{\dot{\phi_{c}}}{\phi_{d}}+B_{cd}\outerproduct*{\phi_{c}}{\dot{\phi_{d}}}\right)-\mathcal{{L}}[\hat{\rho}]\right\\}\ket{\phi_{a}}B_{ab}.\end{split}$
(29)

Inserting the definition of $\dot{\boldsymbol{B}}$ from Eq. (27), we obtain

$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
(30)
$\displaystyle\begin{split}&=\sum\limits_{ab}\bra{\partial_{\alpha_{k}}\phi_{b}}\Big{\\{}\sum\limits_{cd}\big{(}[\boldsymbol{S}^{-1}\boldsymbol{L}\boldsymbol{S}^{-1}\\\
&\qquad\qquad\qquad-\boldsymbol{S}^{-1}\boldsymbol{\tau}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{\tau}^{\dagger}\boldsymbol{S}^{-1}]_{cd}\outerproduct*{\phi_{c}}{\phi_{d}}\\\
&\qquad+B_{cd}\outerproduct*{\dot{\phi_{c}}}{\phi_{d}}+B_{cd}\outerproduct*{\phi_{c}}{\dot{\phi_{d}}}\big{)}-\mathcal{{L}}[\hat{\rho}]\Big{\\}}\ket{\phi_{a}}B_{ab}\end{split}$
(31)

By using the definition of the projector $\hat{P}$ in Eq. (19), we have

$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
(32)
$\displaystyle\begin{split}&=\sum\limits_{ab}\bra{\partial_{\alpha_{k}}\phi_{b}}\Big{\\{}\sum\limits_{cd}\big{(}\hat{\mathds{1}}-\hat{P}\big{)}\ket*{\dot{\phi}_{c}}B_{cd}S_{da}\\\
&\qquad\qquad+\big{(}\hat{P}-\hat{\mathds{1}}\big{)}\mathcal{{L}}[\hat{\rho}]\ket{\phi_{a}}\Big{\\}}B_{ab}\end{split}$
(33)
$\displaystyle\begin{split}&=\sum\limits_{abcl}\bra{\partial_{\alpha_{k}}\phi_{b}}[\hat{\mathds{1}}-\hat{P}]\ket*{\partial_{\alpha_{l}}\phi_{c}}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]_{cb}\dot{\alpha}_{l}\\\
&\qquad\qquad-\bra{\partial_{\alpha_{k}}\phi_{b}}[\hat{\mathds{1}}-\hat{P}]\mathcal{{L}}[\hat{\rho}]\ket{\phi_{a}}B_{ab}\end{split},$
(34)

where in the last step we have inserted the definition of
$\ket*{\dot{\phi}_{c}}$. Importantly, since
$\bra{\partial_{\alpha_{k}}\phi_{m}}[\hat{\mathds{1}}-\hat{P}]\ket*{\partial_{\alpha_{l}}\phi_{n}}=0$
for $k\neq l$, as a variation of $\alpha_{k}$ only affects states in mode $k$,
we can drop the sum over $l$ in Eq. (LABEL:eq:app_alpha_5). Defining

$\displaystyle[\boldsymbol{C_{0}}]_{mn}^{k}$
$\displaystyle=\bra{\partial_{\alpha_{k}}\phi_{m}}[\hat{\mathds{1}}-\hat{P}]\ket*{\partial_{\alpha_{k}}\phi_{n}},$
(35) $\displaystyle[\boldsymbol{Y_{0}}]_{mn}^{k}$
$\displaystyle=\bra{\partial_{\alpha_{k}}\phi_{m}}[\hat{\mathds{1}}-\hat{P}]\mathcal{{L}}[\hat{\rho}]\ket{\phi_{n}},$
(36)

we thus obtain

$\displaystyle\text{Tr}\left\\{{\partialderivative{\hat{\rho}}{\alpha_{k}^{*}}}(\dot{\hat{\rho}}-\mathcal{L}[\hat{\rho}])\right\\}$
(37)
$\displaystyle=\sum\limits_{abc}\boldsymbol{C_{0}}_{bc}^{k}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]_{cb}\dot{\alpha}_{k}-\boldsymbol{Y_{0}}_{ba}^{k}\boldsymbol{B}_{ab}$
(38)
$\displaystyle=\text{Tr}\\{\boldsymbol{C_{0}}^{k}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}\dot{\alpha}_{k}-\text{Tr}\\{\boldsymbol{Y_{0}}^{k}\boldsymbol{B}\\}.$
(39)

Finally, the equation of motion for $\dot{\alpha}_{k}$ reads:

$\dot{\alpha}_{k}=\frac{\text{Tr}\\{\boldsymbol{Y_{0}}^{k}\boldsymbol{B}\\}}{\text{Tr}\\{\boldsymbol{C_{0}}^{k}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}}.$
(40)

We hence have an implicit ordinary differential equation defined by Eqs. (5),
(6), reading

$\displaystyle\boldsymbol{\dot{B}}$
$\displaystyle=\boldsymbol{S}^{-1}\boldsymbol{L}\boldsymbol{S}^{-1}-\boldsymbol{S}^{-1}\boldsymbol{\tau}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{\tau}^{\dagger}\boldsymbol{S}^{-1},$
(41) $\displaystyle\dot{\alpha_{k}}$
$\displaystyle=\frac{\text{Tr}\\{\boldsymbol{Y_{0}}^{(k)}\boldsymbol{B}\\}}{\text{Tr}\\{\boldsymbol{C_{0}}^{(k)}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}}.$
(42)

##  II Operator decomposition

We now detail the decomposition of operators and the construction of the
overlap matrix for the coherent-state ladder basis and the cat-ladder basis.

###  II.1 Coherent-state ladder

For the coherent-state ladder with
$\mathcal{S}=\\{\lVert{\alpha,n}\rangle={{}\hat{a}^{\dagger}}^{n}\lVert{\alpha}\rangle,n=0\dots
N\\}$, we can represent combinations of creation and annihilation operators
as:

$\begin{split}\hat{a}\lVert{\alpha,n}\rangle&=\hat{a}{{}\hat{a}^{\dagger}}^{n}\lVert{\alpha}\rangle=({{}\hat{a}^{\dagger}}^{n}\hat{a}+n{{}\hat{a}^{\dagger}}^{n-1})\lVert{\alpha}\rangle\\\
&=\alpha\lVert{\alpha,n}\rangle+n\lVert{\alpha,n-1}\rangle.\end{split}$ (43)

Hence with the overlap matrix $\boldsymbol{S}$ with elements
$S_{mn}=\langle\alpha,m\rVert\alpha,n\rangle$, we have

$\langle{\alpha,m}\rVert\hat{a}\lVert{\alpha,n}\rangle=\alpha
S_{m,n}+nS_{m,n-1}.$ (44)

Representations of other operators are obtained similarly.

We have the following recursion relation for the overlap matrix
$\boldsymbol{S}$:

$\displaystyle S_{mn}=$
$\displaystyle\langle{\alpha,m-1}\rVert\hat{a}\hat{a}^{\dagger}\lVert{\alpha,n-1}\rangle$
(45) $\displaystyle\begin{split}=&S_{m-1,n-1}+(m-1)(n-1)S_{m-2,n-2}\\\
&+|\alpha|^{2}S_{m-1,n-1}+\alpha(m-1)S_{m-2,n-1}\\\
&+\alpha^{*}(n-1)S_{m-1,n-2}.\end{split}$ (46)

This allows to recursively calculate the matrix elements $S_{mn}$.

###  II.2 Cat ladder

In the non-orthonormal basis of Bargmann cat states, we can decompose the
annihilation operator $\hat{a}$ in the subsystem decomposition by calculating
its action on a basis state:

$\hat{a}\lVert{\mathcal{C}_{\alpha,n}^{\pm}}\rangle=\alpha\lVert{\mathcal{C}_{\alpha,n}^{\mp}}\rangle+n\lVert{\mathcal{C}_{\alpha,n-1}^{\mp}}\rangle,$
(47)

and hence in a modular subsystem decomposition [56] the annihilation operator
reads

$\hat{a}=\hat{X}\otimes(\tilde{b}+\alpha).$ (48)

With the inner product matrix

$S_{mn}^{\mu\nu}=\langle{\mathcal{C}_{\alpha,m}^{\mu}}\lVert{\mathcal{C}_{\alpha,n}^{\nu}}\rangle,$
(49)

where indices $\mu,\;\nu\in\\{+,-\\}$ label the parity of the state, we have

$\langle{\mathcal{C}_{\alpha,m}^{\mu}}\rVert\hat{a}\lVert{\mathcal{C}_{\alpha,n}^{\nu}}\rangle=\alpha
S_{mn}^{\mu,\bar{\nu}}+nS_{m,n-1}^{\mu,\bar{\nu}},$ (50)

where $\bar{\nu}$ indicates a flip of parity index $\nu$. Analogous, we can
obtain expressions for the representation of other operators in this basis.

In the cat-ladder basis, the recursion relation to calculate the inner product
matrix elements $S_{mn}^{\mu\nu}$ is analogous to Eqs. (45),(46):

$\displaystyle S_{mn}^{\mu\nu}=$
$\displaystyle\langle{\mathcal{C}_{\alpha,m-1}^{\bar{\mu}}}\rVert\hat{a}\hat{a}^{\dagger}\lVert{\mathcal{C}_{\alpha,n-1}^{\bar{\nu}}}\rangle$
(51)
$\displaystyle\begin{split}=&S_{m-1,n-1}^{\bar{\mu}\bar{\nu}}+(m-1)(n-1)S_{m-2,n-2}^{\mu\nu}\\\
&+|\alpha|^{2}S_{m-1,n-1}^{\mu\nu}+\alpha(m-1)S_{m-2,n-1}^{\mu\nu}\\\
&+\alpha^{*}(n-1)S_{m-1,n-2}^{\mu\nu}.\end{split}$ (52)

With the overlap matrix for the cat-qubit

$S_{00}^{\mu\nu}=\begin{pmatrix}4\cosh{\left(|\alpha|^{2}\right)}&0\\\
0&4\sinh{\left(|\alpha|^{2}\right)}\end{pmatrix},$ (53)

$S_{mn}^{\mu\nu}$ can be recursively computed for $m,n>0$.

##  III Numerical implementation

We here sketch a numerical implementation of the method presented in this
paper. We first consider a single bosonic mode in the cat-ladder basis. Since
calculating a first-order derivative with respect to $\alpha$ directly
corresponds to the application of $\hat{a}^{\dagger}$ and parity conjugation
in the chosen basis, we can easily relate derivatives to block shifts in the
inner product matrix:

$\langle{\mathcal{C}_{\alpha,m}^{\mu}}\lVert{\partial_{\alpha}\mathcal{C}_{\alpha,n}^{\nu}}\rangle=\langle{\mathcal{C}_{\alpha,m}^{\mu}}\rVert\hat{a}^{\dagger}\lVert{\mathcal{C}_{\alpha,n}^{\bar{\nu}}}\rangle=S_{m,n+1}^{\mu,\nu}.$
(54)

We can then calculate the matrices $\mathbf{C}_{0}$ and $\mathbf{Y}_{0}$,
defined by Eqs. (35), (36) in terms of block-shifted matrices. By denoting
block-shifts of matrices in the form
$[\boldsymbol{A}^{(l)}]_{mn}^{\mu\nu}=[\boldsymbol{A}]_{m+1,n}^{\mu\nu}$,
$[\boldsymbol{A}^{(r)}]_{mn}^{\mu\nu}=[\boldsymbol{A}]_{m,n+1}^{\mu\nu}$, and
$[\boldsymbol{A}^{(l,r)}]_{mn}^{\mu\nu}=[\boldsymbol{A}]_{m+1,n+1}^{\mu\nu}$
for any Matrix $\boldsymbol{A}$ with elements
$[\boldsymbol{A}]_{mn}^{\mu\nu}$, we can rewrite $\boldsymbol{C_{0}}$ and
$\boldsymbol{Y_{0}}$ as

$\displaystyle\mathbf{C}_{0}$
$\displaystyle=\boldsymbol{S}^{(l,r)}-\boldsymbol{S}^{(l)}\boldsymbol{S}^{-1}\boldsymbol{S}^{(r)}$
(55) $\displaystyle\mathbf{Y}_{0}$
$\displaystyle=\sum_{p}\left(\mathbf{A}^{(l)}_{p}\mathbf{BD}_{p}-\mathbf{S}^{(l)}\mathbf{S}^{-1}\mathbf{A}_{p}\mathbf{BD}_{p}\right)$
(56)
$\displaystyle=\mathbf{L}^{(l)}-\mathbf{S}^{(l)}\mathbf{S}^{-1}\mathbf{L}.$
(57)

Here we have expanded the Liouvillian matrix $\boldsymbol{L}$ in its operator-
sum representation,
$\boldsymbol{L}=\sum_{p}\boldsymbol{A}_{p}\boldsymbol{B}\boldsymbol{D}_{p}$,
with left and right Kraus operator matrices $\boldsymbol{A}_{p}$ and
$\boldsymbol{D}_{p}$.

Note that matrix $\mathbf{C}_{0}$, which can be interpreted as a quantum-
geometric tensor, is not full rank, but only has two diagonal non-zero
elements associated with the highest moment of $\hat{a}^{\dagger}$ included in
the basis (one for each parity sector). Therefore, for any initial density
matrix that has no population for the highest moment included in the basis
$\text{Tr}\\{\boldsymbol{C_{0}}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}$
will be zero. As a result $\dot{\alpha}$ will be zero until the density matrix
$\boldsymbol{B}$ populates the highest moment included in the basis ansatz.
Intuitively, this can be understood in the following way: For our time-
dependent variational ansatz, the equations of motion describe a dynamics in
which the density matrix $\boldsymbol{B}$ evolves into a full-rank matrix.
Upon reaching the corner (i.e. the highest moment of $\hat{a}^{\dagger}$
included in the basis), the basis re-centers according to the evolution of
$\alpha$. In the numerical simulations, we therefore evolve $\alpha$ only if
$|\text{Tr}\\{\boldsymbol{C_{0}}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\}|>\epsilon$,
with a small threshold $\epsilon$, otherwise we set $\dot{\alpha}=0$ and
evolve only $\boldsymbol{B}$.

For multiple bosonic modes, the numerical implementation can be
straightforwardly extended. We first note that due to the choice of the basis,
the inner product matrix $\boldsymbol{S}$ can be expressed as a tensor product
of inner product matrices of each bosonic mode, i.e.
$\boldsymbol{S}=\otimes_{k}\boldsymbol{S}_{k}$, also reducing the computation
of the inverse of $\boldsymbol{S}$ to the inverses of $\boldsymbol{S}_{k}$.
Using block-shifts of matrices for mode $k$, we can express left and right
derivatives of matrices with respect to $\alpha_{k}$. Analogous to Eqs. (55),
(57), we obtain

$\displaystyle\mathbf{C}_{0}^{(k)}$
$\displaystyle=\boldsymbol{S}^{(l_{k},r_{k})}-\boldsymbol{S}^{(l_{k})}\boldsymbol{S}^{-1}\boldsymbol{S}^{(r_{k})}$
(58) $\displaystyle\mathbf{Y}_{0}^{(k)}$
$\displaystyle=\sum_{p}\left(\mathbf{A}^{(l_{k})}_{p}\mathbf{BD}_{p}-\mathbf{S}^{(l_{k})}\mathbf{S}^{-1}\mathbf{A}_{p}\mathbf{BD}_{p}\right)$
(59)
$\displaystyle=\mathbf{L}^{(l_{k})}-\mathbf{S}^{(l_{k})}\mathbf{S}^{-1}\mathbf{L},$
(60)

where here the indices $l_{k}$ and $r_{k}$ indicate block-shifts associated to
mode $k$.

The numerical implementation for the coherent-state ladder basis can be
straightforwardly obtained by removing the parity indices $\mu$, $\nu$.

##  IV Parametrizing the basis using higher-order processes

In addition to the parametrization of the basis by coherent-state amplitude
$\alpha$ alone, we can include higher-order processes, such as squeezing and
other higher-order processes. For a single bosonic mode, we define the
following un-normalized state including up to $K$-photon moments,

$\lVert{\alpha^{(1)},\dots,\alpha^{(K)}}\rangle\equiv\prod\limits_{n=1}^{K}e^{\alpha^{(n)}{{}\hat{a}^{\dagger}}^{n}}\ket{0}.$
(61)

We can naturally relate derivatives of the parameters $\alpha^{(n)}$ to the
application of ${{}\hat{a}^{\dagger}}^{n}$:

$\partialderivative{\alpha^{(n)}}\lVert{\alpha^{(1)},\dots,\alpha^{(K)}}\rangle={{}\hat{a}^{\dagger}}^{n}\lVert{\alpha^{(1)},\dots,\alpha^{(K)}}\rangle.$
(62)

For a single bosonic mode, the equation of motion for $\dot{\alpha}^{(n)}$
then becomes

$\dot{\alpha}^{(n)}=[\boldsymbol{C^{-1}}\boldsymbol{Y}]_{n},$ (63)

with

$\displaystyle[\boldsymbol{C}]_{kl}$
$\displaystyle=\text{Tr}\\{\boldsymbol{C_{0}}^{(k,l)}[\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}]\\},$
(64) $\displaystyle\boldsymbol{Y}_{k}$
$\displaystyle=\text{Tr}\\{\boldsymbol{Y_{0}}^{(k)}\boldsymbol{B}\\},$ (65)
$\displaystyle[\boldsymbol{C_{0}}^{(k,l)}]_{mn}$
$\displaystyle=\matrixelement{\partialderivative{\phi_{m}}{\alpha^{(k)}}}{[\hat{\mathds{1}}-\hat{P}]}{\partialderivative{\phi_{n}}{\alpha^{(l)}}}$
(66) $\displaystyle[\boldsymbol{Y_{0}}^{(k)}]_{mn}$
$\displaystyle=\matrixelement{\partialderivative{\phi_{m}}{\alpha^{(k)}}}{[\hat{\mathds{1}}-\hat{P}]\mathcal{{L}}[\hat{\rho}]}{\phi_{n}}.$
(67)