# Limited Associativity Makes Concurrent Software Caches a Breeze

Dolev Adas
Computer Science
Technion
<EMAIL_ADDRESS>Gil Einziger
Computer Science
Ben-Gurion University of the Negev
<EMAIL_ADDRESS>Roy Friedman
Computer Science
Technion
<EMAIL_ADDRESS>
###### Abstract

Software caches optimize the performance of diverse storage systems, databases
and other software systems. Existing works on software caches automatically
resort to fully associative cache designs. Our work shows that limited
associativity caches are a promising direction for concurrent software caches.
Specifically, we demonstrate that limited associativity enables simple yet
efficient realizations of multiple cache management schemes that can be
trivially parallelized. We show that the obtained hit ratio is usually similar
to fully associative caches of the same management policy, but the throughput
is improved by up to x$5$ compared to production-grade caching libraries,
especially in multi-threaded executions.

##  1 Introduction

_Caching_ is a fundamental design pattern for boosting systems’ performance.
Caches store part of the data in a fast and close memory to the program’s
execution. Accessing data from such a memory, called a _cache_ , reduces
waiting times, thereby improving running times. A _cache hit_ means that an
accessed data item is already located in the cache; otherwise, it is a _cache
miss_. Usually, only a small fraction of the data can fit inside a given
cache. Hence, a _cache management scheme_ decides which data items should be
placed in the cache to maximize the _cache hit ratio_ , i.e., the ratio of
cache hits to all accesses. More so, such a management scheme should be
lightweight; otherwise, its costs would outweigh its benefits.

At coarse granularity, we can distinguish between _hardware caches_ and
_software caches_. A CPU’s hardware cache uses fast SRAM memory as a cache for
the slower main memory, which is made of DRAM [22]. In addition, the main
memory (DRAM) is often utilized as a cache for secondary storage (disks and
SSDs) or network/remote storage. Similarly, most storage systems include a
software-managed cache, e.g., the Velocix media platform stores video chunks
in a multi-layer cache [8]. Additional popular examples include storage
systems like Redis [7], Cassandra [2], HBase [3] and Accumulo [1], graph
databases like DGraph [5] and neo4j [6], etc. Web caching is another popular
example.

An important aspect of the cache management scheme is to select a _victim_ for
eviction whenever there is not enough room in the cache to admit a newly
accessed data item. In _fully associative_ caches, the cache management
schemes can evict any of the cached items, as illustrated in Figure 1. In
contrast, _limited associativity_ caches restrict the selection of data items
that may be evicted. Specifically, limited associativity caches are portioned
into independent _sets_. Each set consists of $k$ _ways_ (places), each of
which can store a data item. When admitting an item to a limited associativity
cache, one must select the cache victim from the same set as the newly
admitted item. In such caches, known as _k-way set associative caches_ ,
items’ IDs are mapped into sets using a hash function and each set is an
_independent_ sub-cache of size $k$.

In hardware, limited associativity caches are simpler to implement, consume
less power to operate and are often significantly faster than fully
associative designs. Thus, the vast majority of hardware-based caches have
limited associativity. In contrast, non-distributed software-based caches are
almost always fully associative. Specifically, cache implementations usually
employ hash tables as building blocks. This is often so obvious that papers
describing software caches do not even explicitly mention that their
algorithms are fully associative [14, 17, 33, 41]. A possible reason for this
lack of previous consideration stems from the understanding that fully
associative caches should yield higher hit ratios than the same sized limited
associativity caches. Further, it is relatively easy to implement such caches
from fully associative building blocks such as concurrent hash tables.

(a) Fully Associative Cache

(b) Limited Associativity Cache

Figure 1: A schematic overview of a fully associative cache (a) and a limited
associativity cache (b). In a fully associative cache, Alice and Bob operate
on a shared resource and require some form of synchronization between them. In
the limited associativity case, the cache is partitioned into many independent
sub-caches, also known as sets, and updates are assigned to sets through
hashing. Most times, Alice and Bob operate on different sets requiring no
synchronization at all.

However, our work argues that software caches significantly benefit from a
limited associativity design due to the following advantages: First, limited
associativity caches are embarrassingly parallel and allow for concurrent
operations on different sets without any form of synchronization (as
illustrated in Figure 1). In contrast, it is far from trivial to introduce
concurrency to most fully associative cache policies as often these cache
management policies rank cached data items and then choose the “worst” item
for eviction, c.f., [28, 39]. Thus, parallel actions may select the same cache
victim or update the ranks of cached data items. For example, in LRU, we get
contention on the head of the LRU list [20]. Aside from parallelism, cache
implementations rely on elaborate data structures that often require an
excessive number of pointers. Also, hash tables almost always have a constant
fraction of unused space. In comparison, limited associativity caches offer a
denser representation of data. When $k$ is reasonably small, we can scan all
the items in the set without requiring any auxiliary data structures.

The above limitations of full associativity have motivated the development of
reduced accuracy management policies such as Clock [15, 25], sampled LRU [7]
and sampled dynamic priority based policies [14, 12]. As elaborated in this
work, we claim that limited associativity is a beneficial design alternative
compared to Clock and sampling, especially for parallel systems.

###  1.1 Contributions

In this work, we explore and promote limited associativity designs for
software caches. Specifically, our first contribution is the implementation of
several leading cache replacement policies in a limited associativity manner.
As we show, limited associativity leads to simplified data structures and
algorithms. These can be trivially implemented in static memory with minimal
overhead and operate in constant time. Further, these realizations trivially
support wait-free concurrent manipulation and avoid hot-spots.

Our second contribution is an evaluation of the hit-ratio obtained by multiple
cache management policies on real traces for $k$-way set associative caches
and for fully associative caches. Our findings indicate that, even for
relatively small values of $k$, the difference between the hit-ratio obtained
by a $k$-way associative scheme and the corresponding fully associative scheme
is marginal.

Since the runtime of each operation depends on the size of each set, low
associativity is preferred for speed. Low associativity also implies more
independent sets which reduces the contention on each set, and allows for
better parallelism. However, this is a tradeoff as associativity also impacts
the hit-ratio and if we set the associativity too low then we would suffer a
non-negligible drop in hit-ratio. Our work demonstrates that a reasonable
associativity value of $8$ provides the best of both worlds as it yields a
very similar hit-ratio to fully associative caches, while being considerably
faster.

These results echo previous studies from hardware caching [22, 24], where
limited associativity was found to be nearly as effective as full
associativity. However, those studies were mostly conducted using hardware
workloads, which tend to exhibit higher locality than software workloads.

We stress that limited associativity is a design principle rather than a cache
management policy. Hence, we can adjust most existing policies to their
respective limited associativity version with small changes. This assertion is
especially straightforward for sampling based policies such as sampled LRU
[7], Hyperbolic [14] and LHD [12] that randomly select a small number of
cached items and choose the victim only from these items. Notice that limited
associativity only requires computing a single hash function per cache miss
compared to invoking the PRNG multiple times in sampled approaches. The
sampled approach also requires accessing multiple random memory locations,
which is not friendly for the hardware cache. In contrast, in limited
associativity, we access a short continuous region of memory.

Further, contemporary cache management schemes, including ARC [33], LIRS [26],
FRD [38] and W-TinyLFU [17, 16] maintain two or more cache regions, each of
which handled in a fully associative manner. We argue that each cache region
could be treated as a corresponding limited associativity region for these
schemes.

Last, we evaluate the speedups obtained by our concurrent implementations
compared to leading alternatives. These demonstrate the performance benefit of
limited associativity designs for software caches in modern multi-core
computers.

##### Paper Organization

Section 2 surveys related work and provides the background necessary to
position our work within the broader context of relevant designs. Section 3
suggests several ways to implement $k$-way set associative caches efficiently
in software. Section 5 shows an evaluation of the hit ratio obtained by
limited associativity caches for a wide selection of caching algorithms and of
their throughput compared to the Guava and Caffeine Java caching libraries. We
conclude with a discussion in Section 6.

##  2 Related Work

###  2.1 Cache Management Policies

A cache management scheme’s design often involves certain (possibly implicit)
assumptions on the characteristics of the “typical” workload. For example, the
_Least Recently Used (LRU)_ policy always admits new items to the cache and
evicts the least recently used items [22]. LRU works well for _recency biased
workloads_ , where recently accessed items are most likely to be accessed
again.

Alternatively, the _Least Frequently Used (LFU)_ policy, also called _Perfect
LFU_ , assumes that the access distribution is fixed over time, meaning that
frequently accessed data items are more likely to be reaccessed. Hence,
Perfect LFU evicts the least frequently used item from the cache and admits a
new one if it is more frequent than the cache victim. For synthetic workloads
with static access distributions, Perfect LFU is the optimal cache management
policy [27]. Real workloads are dynamic, so practical LFU policies often apply
aging mechanisms to the items’ frequency. Such mechanisms can be calculating
frequency on sliding windows [27], using exponential decay [10, 9], or
periodic aging [17].

Operation complexity is another essential factor in cache design. Both LRU and
LFU need to maintain a priority queue or a heap, and can be implemented in
constant time [29, 34, 13]. LRU and variants of LFU also maintain a hash table
for fast access to items and their metadata. Another limitation of priority
queue based LRU implementations is the high contention on the head of queue
due to the need to update it on each cache access (both hits and misses).

The realization that a good cache management mechanism needs to combine
recency and frequency led to the development of more sophisticated approaches
[20]. LRU-K [37] maintains the last $K$ occurrences of each item and orders
all items based on the recency of their last $K^{\mathrm{th}}$ access. The
drawbacks are significant memory and computational overheads. Also, most
benefit comes when $K=2$, with quickly diminishing returns for larger values
of $K$.

_Low Inter-reference Recency Set_ (LIRS) [26] is a page replacement algorithm
that attempts to directly predict the next occurrence of an item using a
metric named _reuse distance_. To realize this, LIRS maintains two cache
regions and a large number of ghost entries. FRD [38] can be viewed as a
practical variant of LIRS, which overcomes many of the implementation related
limitations of LIRS.

_Window-TinyLFU_ (W-TinyLFU) [17] also balances between recency and frequnecy
using two cache regions, a _window cache_ and the _main cache_. Yet, it
controls which items are admitted to the main cache using an optimized
counting Bloom filter [19] based admission filter and thereby avoids the need
to maintain ghost entries. The relative size of the two caches is dynamically
adjusted at runtime [16].

As we mentioned before, what is common to ARC, CAR, LIRS, SLRU, 2Q, FRD and
W-TinyLFU is that they all maintain two or more caches regions, each of which
is maintained as a fully associative cache. We promote implementing each such
cache region as a limited associativity cache.

Hyperbolic caching [14] logically maintains a dynamically evolving priority
for each cached item, which is the number of times that item has been accessed
divided by the duration of time since the item was inserted into the cache. To
make this practical and efficient, priorities are only computed during the
eviction and only for a small sample of items. Among the sampled items, the
one whose priority is smallest becomes the victim. As we discuss later, this
policy can be trivially implemented using an associative design.

###  2.2 Associativity in Caches

The only limited associativity software cache we are aware of is HashCache
[11]. Yet, that work focused on very memory constraint Web caches in which the
meta-data cannot fit in main memory and the goal is to minimize disk accesses
to the meta-data. Further, it only considered the LRU replacement policy and
only very particular workloads. In this paper, we perform a more systematic
study of limited associativity’s effectiveness to software cache design.

Current caching libraries use fully associative caches. Yet, in principle, any
distributed cache (e.g., memcacheD) can be viewed as relying on limited
associativity, since each node manages its caches independently to the others.
When an item is added to a given cache, the victim must be selected from the
same host. The main difference though is that in distributed caching, each
host is typically responsible for a very large range of objects (equivalent to
having a very large $k$ value). In contrast, in $k$-way associativity, our
goal, as explained before, is to have a large number of small sets (low $k$
values). Further, the CPU’s L1, L2 and L3 caches are almost always implemented
through limited associativity, with the exception of TCAM based caches. In
hardware, it is clear that fully associative caches require fewer circuits to
implement and that they work faster.

###  2.3 On Hash Tables and Caches

As mentioned above, most fully associative software caches use fully
associative hash tables to quickly find the cached items. In principle, hash
tables incur non-negligible overheads. Specifically, open addressing hash
tables such as the classic Linear Probing, or the modern Hopscotch hash table
[23], are implemented on top of partially full arrays. In all open address
hash tables, a constant fraction of these arrays must remain empty. Thus,
relying on such tables results in underutilized space that could potentially
store useful items. Closed address hash tables such as Java’s concurrent hash
map [4] are another alternative. Yet, they incur significant overheads for
pointers, which results in a suboptimal data density.

In contrast, limited associativity caches can be implemented using the same
arrays as open address hash tables, but we can use 100% of the array’s cells.
In principle, we do not even have to use pointers and achieve a very dense
representation of data. Such a trade-off is potentially attractive since our
caches can store more items within a given memory unit and having more items
increases the cache hit ratio. On the other hand, there is less freedom in
selecting the cache victim, which may reduce the hit ratio.

Let us emphasize that hashtables and caches are not the same. While a
hashtable can be used to store cached objects, one would still need additional
data-structures and algorithms to implement the cache management policy, and
this is the main source of pain in software caches. On the other hand, caches
are allowed to evict items while hashtables retain all the items that were not
explicitly removed.

###  2.4 Parallel Caches

Caches need to operate fast and thus parallelism is a vital consideration.
Specifically, LRU is known to yield high contention on the head of the list as
a simple implementation moves elements to the head of the list on each access.
This contention motivated alternative designs such as Clock [15]. In Clock,
each entry receives a single bit, indicating whether it was recently accessed.
Items are always admitted by Clock and the cache victim is the first item with
an unset Clock bit. The Clock hand clears set bits as it searches for the
cache victim. Although the worst-case complexity of Clock is linear, it
parallelizes well and creates little contention and is therefore used in
practice as an operating system page replacement cache. Another interesting
approach is taken by [18], which uses Bloom filter based approximations of LRU
to enable better concurrency.

Fully associative cache management often relies on a hash table that needs to
be parallelized. Such tables are a shared resource and are usually not
embarrassingly parallel [35], [30]. On the other hand, limited associativity
caches are embarrassingly parallel as actions on different sets are completely
independent and require no synchronization. This property holds regardless of
the (limited associativity) eviction policy.

##  3 Limited Associativity Architecture

A K-Way cache supports two operations: get/read and put/write, as customary in
other caches such as Caffeine [31] and Guava [21]. A get/read operation
retrieves an item’s value from the cache or returns null if this item is not
in the cache. The put/write operation inserts an item into the cache; if the
item already exists in the cache, this operation overwrites its value. Both
types of operations update the item’s metadata that is used by the management
policy.

Figure 2: An example of a 48 items cache organized in a 6-way set associative
design. The cache is broken into 8 independent sets such that each contains a
up to 6 items.

Figure 2 illustrates such cache organization consisting of $8$ independent
sets that are held in an array. Each set is an array of $6$ ways (nodes) that
store items with metadata that may be required by various eviction policies.

Upon insertion of an item to the cache, we use a hash function to determine
its set. If the set is full, we also need to select a victim from the same set
to be evicted. K-Way supports efficient implementations of numerous eviction
policies without relying on expensive auxiliary data structures. For example,
to implement LRU or LFU, we only need a short counter for each cached item in
the set. We select the victim by reading all $K$ items in the set and evicting
the one with the smallest counter. The difference between LRU and LFU is
merely in the way the counter is updated. In LRU, we use it to record the
(logical) timestamp of the last access, whereas in LFU, we use it to count how
many times the entry was accessed in the past. To read an item, we scan the
appropriate set. Once found, we update the corresponding item’s metadata
(updating the access time and frequency) atomically and returning the value.

Clearly, the K-Way time complexity of both get/read and put/write is O($K$),
which is only dependent on the number of ways ($K$). As $K$ is a small
constant number, we treat it as constant time.

Similarly, for Hyperbolic caching [14] we maintain two short counters for each
item, one recording the time it was inserted ($t_{0}$) and the other ($n$)
counting the number of times it was accessed (initialized to $1$). For an
eviction at time $t_{e}$, we scan the $K$ items at the corresponding set and
select the item whose $n/(t_{e}-t_{0})$ value is minimal.

We implemented K-Way in Java. It supports five eviction policies: LRU, LFU,
FIFO, Random and Hyperbolic. The eviction policy is chosen at the constructor.
We have realized several concurrency control variants as detailed below.

#####  $K$ Way Cache Wait Free Array - KW-WFA

Each set is an array of Node references. Algorithm 1 depicts the internal
classes and variables. Pseudo-code for the get operation is given in Algorithm
2 and for the put operation in Algorithm 3. To replace a victim node with a
newly arriving one, we perform a CAS (compare and swap) operation on the
address of this node.

Algorithm 1 Internal KW-WFA classes and fields

1:class Node $\\{$

2:K key;

3:V value;

4:AtomicInteger counter;

5:int index; $\\}$

6:class set $\\{$

7:AtomicReferenceArray$<$Node$<$K, V$>$$>$ setArray;

8:AtomicLong time; $\\}$ // time is only present in LRU

Algorithm 2 KW-WFA get(read) operation

1:function get(K key)

2: int set =(int) hash(key) & (numberOfSets-1);

3: for int i = 0; i $<$ ways; ++i do

4: Node$<$K, V$>$ n = cache[set].setArray.get(i);

5: if n != null && n.key.equals( key) then

6: update(n.counter);

7: return n.value;

8: return null;

Algorithm 3 KW-WFA put(write) operation

1:function put(K key , V value)

2: int set =(int) hash(key) & (numberOfSets-1));

3: for int i = 0; i $<$ ways; ++i do

4: Node$<$K, V$>$ n = cache[set].setArray.get(i);

5: if n != null && n.key.equals( key) then

6: cache[set].setArray.compareAndSet(i,n,update);

7: return ;

8: Node$<$K, V$>$ victim =policy.select(cache[set].setArray, key)

9: if victim != null then

10: Node$<$K,V$>$ updateNode = new
Node$<$K,V$>$(key,value,victim.index,cache[set].readTime());

11: cache[set].setArray.compareAndSet(i,victim,update);

12: else

13: for int i = 0; i $<$ ways; ++i do

14: if Node$<$K, V$>$) cache[set].setArray.get(i)== null then

15: Node$<$K, V$>$ updateNode = new Node$<$K,
V$>$(key,value,victim.index,cache[set].readTime());

16: cache[set].setArray.compareAndSet(i ,null,updateNode);

17: return;

To exchange an entire node atomically, we maintain each set of K-Way as a
reference array. Alas, this has the following shortcoming. To search an item
and select a victim, we must scan the entire corresponding set reading the
metadata of each node. With a reference array implementation, this means
reading multiple different memory locations, one for each item in the set,
which is expensive.

To improve this and enable scans to access continuous memory regions, we offer
a second implementation in which we separated the counters and fingerprints
from the nodes.

##### K Way Cache Wait Free Separate Counters - KW-WFSC

We store an array of counters and an array of fingerprints with the reference
array of the nodes for each set. Algorithm 4 depicts the internal classes and
variables. To insert an item, we first scan the fingerprints array. If a
certain entry’s fingerprint matches the key’s fingerprint, we check the key
inside the corresponding node. If they are the same, we update this node with
the new value and return. If the item is new, we scan the counter array to
find the lowest/highest counter, depending on the eviction policy. We then
replace the victim without accessing the node at all. Pseudo-code for the get
operation is listed in Algorithm 5 and for the put operation can be found in
Algorithm 6.

Algorithm 4 Internal KW-WFSC classes and fields

1:class Node $\\{$

2:K key;

3:V value;

4:AtomicInteger counter;

5:int index; $\\}$

6:class set $\\{$

7:AtomicReferenceArray$<$Node$<$K, V$>$$>$ setArray;

8:AtomicIntegerArray counters;

9:AtomicLongArray fingerPrint; $\\}$

Algorithm 5 KW-WFSC get operation

1:function get(K key)

2: int set =(int) hash(key) & (numberOfSets-1);

3: for int i = 0; i $<$ ways; ++i do

4: Long fp = cache[set].fingerPrint.get(i);

5: if fp.equals((Long)key) then

6: Node$<$K, V$>$ n = cache[set].setArray.get(i);

7: if n != null && n.key.equals( key) then

8: update(n.counter);

9: return n.value;

10: return null;

Algorithm 6 KW-WFSC put operation

1:function put($K\ key,\ V\ value$)

2: int set =(int) hash(key) & (numberOfSets-1));

3: for int i = 0; i $<$ ways; ++i do

4: Long fp = cache[set].fingerPrint.get(i);

5: if fp.equals((Long)key) then

6: Node$<$K, V$>$ n = cache[set].setArray.get(i);

7: if n != null && n.key.equals( key) then

8: Node$<$K, V$>$ update = new Node$<$K, V$>$(key, value, i);

9: cache[set].setArray.compareAndSet(i, n, update);

10: return;

11: Node$<$K,V$>$ victim =policy.select(cache[set].counters, key);

12: Node$<$K,V$>$ update = new
Node¡K,V¿(key,value,victim.index,cache[set].readTime());

13: if cache[set].setArray.compareAndSet(victim.index ,victim,update); then

14: cache[set].fingerPrint.set(victim,index,(Long)key);

15: cache[set].counters.set(victim,index, 0);

#####  $K$ Way Cache Lock Set - KW-LS

In our third implementation, we utilize one lock per set. Algorithm 7 depicts
the internal classes and variables. To read an item from the set, we lock the
set and scan it. If the item is in the set, we try to change the lock to a
write lock to update the counter, return the value and finish. We treat writes
in a similar manner. Pseudo-code for the get operation is given in Algorithm 8
and for the put operation in Algorithm 9.

Algorithm 7 Internal KW-LS classes and fields

1:class Node $\\{$

2:$\ \ K\ key;$

3:$\ \ V\ value;$

4:$\ \ Integer\ counter;\\}$

5:$\ \ int\ index;\\}$

6:class set $\\{$

7:$\ \ Node<K,V>[]setArray;$

8:$\ \ StampedLock\ lock;\\}$

Algorithm 8 KW-LS get operation

1:function get(K key)

2: int set =(int) hash(key) & (numberOfSets-1));

3: long stamp =cache[set].lock.readLock();

4: for int i = 0; i $<$ ways; ++i do

5: Node$<$K, V$>$ n = cache[set].setArray.get(i);

6: if n != null && n.key.equals( key) then

7: long stampConvert=cache[set].lock.tryConvertToWriteLock(stamp);

8: if stampConvert == 0 then

9: cache[set].lock.unlockRead(stamp);

10: return n.value;

11: update(n.counter);

12: cache[set].lock.unlockWrite(stampConvert);

13: return n.value;

14: cache[set].lock.unlockRead(stamp);

15: return null;

Algorithm 9 KW-LS put operation

1:function put( K key , V value)

2: int set =(int) hash(key) & (numberOfSets-1));

3: long stamp =cache[set].lock.readLock();

4: for int i = 0; i $<$ ways; ++i do

5: Node$<$K, V$>$ n = cache[set].setArray[i];

6: if n != null && n.key.equals( key) then

7: long stampConvert=cache[set].lock.tryConvertToWriteLock(stamp);

8: if stampConvert == 0 then

9: cache[set].lock.unlockRead(stamp);

10: return ;

11: n.value=value;

12: update(n.counter)

13: cache[set].lock.unlockWrite(stampConvert);

14: return;

15: Node$<$K, V$>$ victim =policy.select(cache[set].setArray, key)

16: int victimIndex = 0

17: if victim != null then

18: victimIndex = victim,index

19: else

20: for int i = 0; i $<$ ways; ++i do

21: if Node$<$K, V$>$ cache[set].setArray.get(i)== null then

22: victimIndex =i

23: Node$<$K, V$>$ updateNode = new Node$<$K,
V$>$(key,value,victim.index,cache[set].readTime());

24: cache[set].setArray[victimIndex].key=key;

25: cache[set].setArray[victimIndex]value.=value;

26: cache[set].setArray[victimIndex].counter=0;

27: cache[set].lock.unlockWrite(stampConvert);

28: return;

##  4 Formal Analysis

Assume that a fully associative cache policy would like to store $C$ items in
the cache. We can guarantee that any specific $C$ items can be kept in a
larger limited associativity cache. Specifically, if we denote the number of
sets $n$, then if $C\geq n\log(n)$ we get that the maximum load (among the $C$
items) is: $\frac{C}{n}+\Theta\left({\sqrt{\frac{C\log(n)}{n}}}\right)$.

The next theorem guarantees that a K-way set associative cache can mimic a
(smaller) fully associative cache with high probability. Thus, such caches are
capable of comparable performance as their fully associative counterparts. Our
theorem uses the following Chernoff bound (Theorem 4.4, page 66) from [36]:
Let $X_{1},X_{2},...X_{n}$ be independent Poission trails, with
$E(X_{i})=p_{i}$, let $X=\Sigma^{n}_{i=1}X_{i}$ and let $\mu=E(X)$ then:
$\Pr\left[X\geq(1+\delta)E(X)\right]\leq e^{\frac{-E(X)\delta^{2}}{3}}.$

######  Theorem 4.1.

Consider a K-Way set associative cache of size $C^{\prime}\geq 2C$, then the
probability for the cache to be able to store the desired $C^{\prime}$ items
is at most: $\frac{C^{\prime}}{k}e^{\frac{k}{6}}$.

###### Proof.

When storing $C$ items in a $C^{\prime}$ sized limited associativity cache
partitioned into $\frac{C^{\prime}}{k}$ sets, each set receives
$\frac{Ck}{C^{\prime}}=\frac{k}{(1+\delta)}$ of those $C$ items. We look at an
arbitrary set and denote $X_{i}$ the random variable of the $i^{\mathrm{th}}$
item. That is, $X_{i}=1$ if the $i^{\prime}th$ item is placed on the set in
question and zero otherwise. We denote by $X$, the total number of items
(among the $C$) that are placed in the set. By definition, we get that:
$E(X)=\frac{k}{(1+\delta)}$.

From the Chernoff bound we get: $\Pr\left[X\geq(1+\delta)E(X)\right]\leq
e^{\frac{k\delta^{2}}{(1+\delta)3}}$. That is, the probability for each set to
not be able to store one of the items allocated to it is:
$e^{\frac{k\delta^{2}}{(1+\delta)3}}$. We apply the union bound on all sets
and obtain that the probability that none of the sets receives more items than
it can store is at most:
$\frac{C^{\prime}}{k}e^{\frac{k\delta^{2}}{(1+\delta)3}}$. Substituting
$\delta=1$ concludes the proof. ∎

For example, Theorem 4.1 shows that a 64-way cache of size 200k items can
store any desired 100k items with a probability of over $99\%$. Alternatively,
we can store in a 2M sized 128 way set associative cache any 1M items with a
probability of over $99.999\%$.

In practice, even when limited associativity caches miss out on some desired
items, they store enough other items to yield comparable hit ratios as same
size fully associative caches. This is due to the following two reasons:
First, the above bound is not tight. Second, it is quite reasonable that
slight content divergence between caches would still result in similar hit
ratios. That is, the exact cache content is just a mean to obtain good hit
ratio, not a goal by itself.

##  5 Evaluation

Figure 3: Legend for all graphs.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) LRU +TinyLFU

Figure 4: wiki1190322952.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic

Figure 5: P8.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic

Figure 6: P12.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic+TinyLfu

Figure 7: S1.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic+TinyLfu

Figure 8: S3.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic+TinyLfu

Figure 9: OLTP.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) FIFO+TinyLFU

Figure 10: LIRS - multi2.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) FIFO+TinyLFU

Figure 11: multi3.

(a) LRU

(b) LFU +TinyLFU

(c) Product

(d) Hyperbolic

Figure 12: sprite.

(a) LRU

(b) [LFU +TinyLFU

(c) Product

(d) LFU

Figure 13: WebSearch3.

###  5.1 Workloads and Setup

We utilized the following real world workloads:

wiki1190322952 and wiki1191277217:

two traces from Wikipedia containing 10% of the traffic to Wikipedia during
three months starting in September 2007 [43].

sprite:

From the Sprite network file system, which contains requests to a file server
from client workstations for a two-day period by [26].

multi1:

This trace is obtained by executing two workloads, cs, an interactive C source
program examination tool trace, and cpp, a GNU C compiler pre-processor trace
by [26].

multi2:

This is obtained by executing three workloads, cs, cpp and postgres, a trace
of join queries among four relations in a relational database system from the
UC Berkeley, provided by [26].

multi3:

This is obtained by executing four workloads, cs, cpp, glimpse, a text
information retrieval utility, and postgres provided by [26].

OLTP:

A file system trace of an OLTP server [33].

DS1:

A database trace provided by [33].SMALL

S3:

Search engine trace provided by [33].

S1 and S3:

Search engine traces provided by [33].

P8, P12 and P14:

Windows server disc accesses [33].

F1 and F2:

Traces of transaction processing taken from large financial institution. The
trace is provided by the UMass trace repository [44].

W2 and W3:

Search engine traces provided by [44].

In all figures below, the “$k$ ways” and “sampled” lines are K-Way with
set/sample sizes of $4,8,16,32,64$ & $128$. The “fully associative” line
stands for a linked-list based fully associative implementation. The “product”
figures for each trace present Guava cache by Google [21], Caffeine [31] and
segmented Caffeine [32], which is Caffeine divided into 64 segments to enable
Caffeine to benefit from parallelism at the possible expense of reduced hit-
ratios (plain Caffeine serves updates by a single thread).

In all throughput evaluations, we compare K-Way with $k=8$ to sampled
implementations of corresponding methods with sample size = 8, Guava cache,
Caffeine and segmented Caffeine. All evaluated cache implementations are
written in Java. Our K-Way cache uses xxHash [42] as the hash function that
distributes items to sets.

####  5.1.1 Platform

We ran on the following platforms:

  * •

AMD PowerEdge R7425 server with two AMD EPYC 7551 Processors, each with $32$
2.00GHz/2.55GHz cores that multiplex $2$ hardware threads, so in total, this
system supports $128$ hardware threads. The hardware caches include 32K L1
cache, 512K L2 cache and 8192K L3 cache. It has $128$GB DRAM organized in 8
NUMA nodes, 4 per processor.

  * •

Intel Xeon E5-2667 v4 Processor including $8$ 3.20GHz cores with 2 hardware
threads, so this system supports $16$ hardware threads. It has $128$GB DRAM.
The hardware caches include 32K L1 cache, 256K L2 cache and 25600K L3 cache.

####  5.1.2 Methodology

During throughput evaluation, we employed the following typical cache
behavior: For each element in the trace, we first performed a read operation.
If the element was not in the cache, we initiated a write operation of that
element into the cache.

We started with a warm-up for filling the cache with elements not in the
trace. We did the warm-up first in the main thread that inserted elements up
to the cache size and then in each thread by inserting $size/\\#threads$ non
cached elements that are not in the trace. All threads started the test after
the warm-up simultaneously by waiting on a barrier.

We ran each test for a fixed period, counting the number of operations. The
time was calculated in comparison to the size of the trace and is between 1
and 4 seconds. Each trace was run with a different cache size taking into
account the size of the trace and the hit rate. Each point in the graphs is
the mean taken over 11 runs.

###  5.2 Hit Ratio Evaluation

Figure 3 shows the legend for all the hit ratio and throughput graphs shown in
this work. Figures 4, 6, 7, 8, 9, 10, 11, 12 and 13 depict the results for the
traces described above with different eviction and admission policies. For
brevity, not all traces and combinations are shown. The ones removed
essentially exhibit similar results.

The first graphs for all traces is LRU due to its popularity. The second graph
is LFU eviction with TinyLFU admission, first proposed in [17]. LFU is also a
popular policy, but LFU alone is not a very good management policy as it lacks
aging mechanisms and lacks data about non-cached objects, which are added by
the TinyLFU admission mechanism. The third graph includes the real caching
products: Guava [21], Caffeine [31] and the proof-of-concept segmented
Caffeine (with number of independent segments that match the number of threads
tested) [32]. For brevity, the fourth graph presents different (additional)
policies for different traces. We summarise the results below:

##### LRU eviction policy

Recall that LRU evicts the least recently used item (and admits all items). In
the limited associativity model, we first assign a new item to one of the
small-sized sets based on hashing its key and then evict the least recently
used item of that set. As can be seen from the graphs, for most traces, the
level of associativity has minimal impact on the obtained hit ratio.

##### LFU eviction with LFU admission

The results for this policy (in subfigures b) show that associativity has a
minor impact on the hit ratio. Also, observe that the performance of the
sampled approach is similar to that of limited associativity. In some cases
one is slightly better than the other and in others it is the opposite.

##### Products

In these graphs (subfigures c), we compared Guava cache [21], Caffeine [31]
and segmented Caffeine [32]. In the latter, we used a hash function to map
each element to one of the corresponding Caffeine managed segments. Each such
Caffeine instance was constructed with MAX SIZE/$\\#threads$, to hold MAX SIZE
elements in total. As can be observed, the hit-ratio of segmented Caffeine is
nearly identical to that of Caffeine. Also, notice that Caffeine attains
higher hit-ratio than Guava as it uses the climber policy from [16] compared
to some version of LRU in Guava.

##### Hyperbolic Caching

We also show that limited associativity has little effect on Hyperbolic
caching (figs. 6 and 12) and on Hyperbolic caching with the TinyLFU admission
policy (fig. 8). Notice that for Hyperbolic caching, there are slight
differences between sampling and limited associativity. When combining
Hyperbolic with TinyLFU admission, these differences become barely noticeable.

##### Conclusions:

This motivational study concludes that limited associativity only has a minor
impact on the obtained hit ratio.

###  5.3 Throughput Evaluation on Real Traces

Figure 14: F1 trace with cache size of $2^{11}$ elements and duration of run
of 1 second, run on AMD PowerEdge R7425.  Figure 15: S3 trace with cache size
of $2^{19}$ elements and duration of run of 4 second, run on AMD PowerEdge
R7425.  Figure 16: S1 trace with cache size of $2^{19}$ elements and duration
of run of 4 second, run on AMD PowerEdge R7425.  Figure 17: wiki1190322952
trace with cache size of $2^{11}$ elements and duration of run of 1 second,
run on AMD PowerEdge R7425.  Figure 18: OLTP trace with cache size of $2^{11}$
elements and duration of run of 1 second, run on AMD PowerEdge R7425.  Figure
19: F2 trace with cache size of $2^{11}$ elements and duration of run of 1
second, run on Intel Xeon E5-2667.  Figure 20: W3 trace with cache size of
$2^{19}$ elements and duration of run of 4 second, run on Intel Xeon E5-2667.
Figure 21: multi1 trace with cache size of $2^{11}$ elements and duration of
run of 1 second, run on Intel Xeon E5-2667.  Figure 22: multi2 trace with
cache size of $2^{11}$ elements and duration of run of 1 second, run on Intel
Xeon E5-2667.  Figure 23: multi3 trace with cache size of $2^{11}$ elements
and duration of run of 1 second, run on Intel Xeon E5-2667.  Figure 24: sprite
trace with cache size of $2^{11}$ elements and duration of run of 1 second,
run on Intel Xeon E5-2667.  Figure 25: P12 trace with cache size of $2^{17}$
elements and duration of run of 2 second, run on Intel Xeon E5-2667.  Figure
26: wiki1191277217 trace with cache size of $2^{11}$ elements and duration of
run of 1 second, run on Intel Xeon E5-2667.

The throughput measurements, in Millions of Get/Put operations per second,
with different number of threads are shown in Figures 14, 15, 16, 17 and 18
for the AMD PowerEdge R7425 server and figs. 19, 20, 21, 22, 23, 24, 25 and 26
for the Intel Xeon E5-2667 server.

As can be observed, the K-Way variations KW-WFA, KW-WFSC and KW-LS gain the
highest throughput. Caffeine and sampled exhibited the worst worst-case
performance, and it does not improve with additional threads for all traces
except Sprite (fig. 24), which will be addressed shortly.

For Caffeine, this is explained as put operations in Caffeine are performed in
the background by a single thread. Guava is considerably faster than Caffeine
in traces with significant number of misses because it performs put operations
in the foreground in parallel. As for the sampled approach, on every miss it
needs to calculate $K$ random numbers and access these $K$ random locations to
create the sample. Even when $K$ is as small as 8, this already incurs a
significant overhead.

Our limited associativity design therefore serves misses much faster than
sampled (and Caffeine). However, sampled may be quicker in serving a hit since
it only accesses the (single) accessed object’s meta-data, while our limited
associativity design always accesses the entire set. Further, the larger the
cache is, the higher is the number of sets (since $K$ is independent of the
cache size). Hence, the hash function of the limited associativity approach is
likely to yield better load balancing for larger caches.

For Sprite (fig. 24), the hit ratio is high and the caches are small. Hence,
our limited associativity under-performs compared to sampled. Yet, it still
delivers tens of millions of operations per second.

Figure 27: Synthetic trace of read and write, 100% miss ratio, with cache size
of $2^{21}$ elements, and duration of run of 5 second, run on AMD PowerEdge
R7425.

###  5.4 Synthetic Throughput Evaluation

##### 100% miss ratio:

The worst-case performance for any cache is the case of 100% misses. In this
evaluation, we perform a single get and a single put operation for each item
and each item is requested only once. This extreme case helps us characterize
the operational framework of various cache libraries.

Figure 27 shows the attained throughput (measured in Millions of Get/Put
operations per second). As can be observed, Caffeine has the worst worst-case
performance and it does not improve with additional threads. This is explained
as put operations in Caffeine are performed in the background by a single
thread. Guava is considerably faster than Caffeine, because it performs put
operations in the foreground, and next are our own three implementations of
limited associativity caches. Here, the KW-WFSC is the fastest of our
algorithms because it is optimized towards fast replacement of the cache
victim.

##### 100% hit ratio:

Next, we repeat the last experiment but only perform Get operations to mimic
having 100% hit ratio. Our results are in Figure 28. Notice that here,
Caffeine is considerably faster than all the alternatives since its read
operations are simple hash table read operations. Next is Guava, which is very
fast, and our algorithms are last. Interestingly, our algorithms attain very
similar performance to the 100% miss case, as we always scan the set - a full
scan on a miss and half the set on average for a hit. This implies that their
performance is less sensitive to the workload characteristics. The performance
of Caffeine and Guava varies greatly according to the cache effectiveness. We
conclude that the break-even point between our approach vs. Caffeine and Guava
depends on the hit ratio of the trace. So our next evaluation aims to
determine at what hit ratios limited associativity caches are expected to
outperform these established libraries.

Figure 28: Synthetic trace to mimic 100% hit ratio with cache size of $2^{21}$
elements, and duration of run of 5 second, run on AMD PowerEdge R7425.

##### 95% hit ratio:

We continue and simulate 95% hit ratio by performing one put operations for
every 20 read operations. The results in Figure 29 show that in 95% hit ratio
Guava is the fastest algorithm, while Caffeine is fastest for a small number
of threads. Here, again Caffeine does not scale with the number of threads
because the bottleneck still becomes the write buffer.

Figure 29: Synthetic trace to mimic 95% hit ratio with cache size of $2^{21}$
elements, and duration of run of 5 second, run on AMD PowerEdge R7425

##### 90% hit ratio:

Next, we increase the number of writes to 1 in 10 to simulate 90% hit ratio.
The results in Figure 30 show that KW-WFSC is already faster than Guava and
Caffeine when there are enough threads. That is, we conclude that our approach
attains higher throughput than Caffeine and Guava as long as the hit ratio is
less than 90%, which implies our attractiveness for a very large variety of
workloads and cache sizes.

Figure 30: Synthetic trace to mimic 90% hit ratio with cache size of $2^{21}$
elements, and duration of run of 5 second, run on AMD PowerEdge R7425

##### Conclusions:

In our limited associativity design, get and put have similar performance. In
contrast, in sampled, Caffeine and Guava gets are fast and puts are slow. Our
performance is expected to be better than sampled, Caffeine and Guava when the
hit ratio is below 90%, which captures most real traces and cache sizes. Our
approach scales better with the number of threads in almost all scenarios.

##  6 Discussion

Our work argues that software cache associativity is an important design
choice that should be given consideration. We note that limited associativity
caches are easier and simpler to implement than fully associative caches. They
are embarrassingly parallel and have lower memory overheads compared to fully
associative caches. On the other hand, limited associativity may negatively
impact the hit ratio, which is an important quality measure for every cache.

Sampled techniques, such as Sampled LRU of Redis [7] and the sampling
performed by Hyperbolic caching [14], reduce accuracy in favor of speed. This
is because one can always improve the hit ratio by increasing the cache size,
but there are no simple techniques to improve the operation speed of the
cache. On the other hand, caches’ utility is due to performance gaps between
different ways to fetch data. For example, main memory is faster than
secondary storage. Therefore, the potential benefits from caching are based on
the gap between the operation time of the cache and the operation time of the
alternative method. Thus, improving the operation time of caches makes every
cache hit more beneficial to the application (as long as we do not
significantly reduce the hit ratio).

Our work evaluated limited associativity caches on diverse (real) application
workloads over multiple popular cache management policies. We showed that
limited associativity caches attain a similar (and often nearly identical) hit
ratio compared to fully associative caches for many workloads and policies.
Our work also suggested three different implementations of limited
associativity caches. We showed an improvement of up to x$5$ compared to the
Caffeine and Guava libraries that together dominate the Java ecosystem,
especially in multi-threaded settings. Given these findings, we believe that
such caching libraries would benefit from adopting limited associativity cache
designs.

Among our different implementations, when the trace is uniformly separated
without heavy hitters, so the insertions are usually to different buckets, it
is better to use the KW-LS. If there are many reads and not many writes, it is
better to use the KW-WFSC as the reads are continuous in memory and hence
faster than KW-WFA where the array is an array of pointers so the performance
of the reads is compromised. Otherwise, it is best to use the KW-WFA as it
uses only one atomic operation in contrast to KW-WFC with three atomic
operations, which slow down the cache updates.

Interestingly, for hardware caches, the seminal work of Qureshi, Thompson and
Patt [40] has found that to fully optimize the cache hit ratio, one should
vary the associativity of a cache on a per-set basis in response to the
demands of the program. Supporting such dynamic associativity level in
software is straight-forward and could help improve the hit-ratio gaps between
$8$ ways and fully associative in traces like WS2. Exploring this direction is
left for future work.

##### Acknowledgments:

This work was partially funded by ISF grant #1505/16 and the Technion-HPI
research school as well as the Cyber Security Research Center and the Data
Science Research Center at Ben-Gurion University

## References

  * [1] Apache Accumulo.  2020\.
  * [2] Apache Cassandra.  2020\.
  * [3] Apache HBase.  2020\.
  * [4] Java ConcurrentMap.  https://en.wikipedia.org/wiki/Java_ConcurrentMap, 2020.
  * [5] The DGraph Website.  2020\.
  * [6] The Neo4j Website.  2020\.
  * [7] Using Redis as an LRU Cache.  2020\.
  * [8] S. Akhtar, A. Beck, and I. Rimac.  HiFi: A Hierarchical Filtering Algorithm for Caching of Online Video.  In ACM MM, MM, pages 421–430, 2015.
  * [9] M. Arlitt, L. Cherkasova, J. Dilley, R. Friedrich, and T. Jin.  Evaluating Content Management Techniques for Web Proxy Caches.  In In Proc. of the 2nd Workshop on Internet Server Performance, 1999\.
  * [10] M. Arlitt, R. Friedrich, and T. Jin.  Performance Evaluation of Web Proxy Cache Replacement Policies.  Perform. Eval., 39(1-4):149–164, Feb. 2000.
  * [11] A. Badam, K. Park, V. S. Pai, and L. L. Peterson.  HashCache: Cache Storage for the Next Billion.  In NSDI, pages 123–136, 2009.
  * [12] N. Beckmann, H. Chen, and A. Cidon.  LHD: Improving Cache Hit Rate by Maximizing Hit Density.  In NSDI, pages 389–403, Apr. 2018.
  * [13] R. Ben-Basat, G. Einziger, R. Friedman, and Y. Kassner.  Heavy Hitters in Streams and Sliding Windows.  In IEEE INFOCOM, 2016.
  * [14] A. Blankstein, S. Sen, and M. J. Freedman.  Hyperbolic Caching: Flexible Caching for Web Applications.  In ATC, pages 499–511, 2017.
  * [15] F. J. Corbato.  A Paging Experiment with the Multics System.  Technical Report Project MAC Report MAC-M-384, MIT, May 1968.
  * [16] G. Einziger, O. Eytan, R. Friedman, and B. Manes.  Adaptive Software Cache Management.  In Middleware, pages 94–106, 2018.
  * [17] G. Einziger, R. Friedman, and B. Manes.  TinyLFU: A Highly Efficient Cache Admission Policy.  ACM ToS, 2017.
  * [18] B. Fan, D. G. Andersen, and M. Kaminsky.  MemC3: Compact and Concurrent MemCache with Dumber Caching and Smarter Hashing.  In NSDI, pages 371–384, 2013.
  * [19] L. Fan, P. Cao, J. Almeida, and A. Z. Broder.  Summary cache: a scalable wide-area web cache sharing protocol.  IEEE/ACM Trans. Netw., 8(3):281–293, June 2000.
  * [20] R. Friedman and B. Manes.  Tutorial: Designing Modern Software Caches.  FAST, 2020.
  * [21] Google.  Guava: Google Core Libraries for Java.  https://github.com/google/guava, 2016.
  * [22] J. L. Hennessy and D. A. Patterson.  Computer Architecture - A Quantitative Approach (5. ed.).  Morgan Kaufmann, 2012.
  * [23] M. Herlihy, N. Shavit, and M. Tzafrir.  Hopscotch Hashing.  DISC, pages 350–364, 2008.
  * [24] M. D. Hill and A. J. Smith.  Evaluating Associativity in CPU Caches.  IEEE Transactions on Computers, 38(12):1612–1630, Dec 1989.
  * [25] S. Jiang, F. Chen, and X. Zhang.  CLOCK-Pro: An Effective Improvement of the CLOCK Replacement.  In USENIX Annual Technical Conference, ATEC, pages 35–35, 2005\.
  * [26] S. Jiang and X. Zhang.  LIRS: An Efficient Low Inter-reference Recency Set Replacement Policy to Improve Buffer Cache Performance.  ACM SIGMETRICS, pages 31–42, 2002.
  * [27] G. Karakostas and D. N. Serpanos.  Exploitation of Different Types of Locality for Web Caches.  In ISCC, 2002.
  * [28] R. Karedla, J. S. Love, and B. G. Wherry.  Caching Strategies to Improve Disk System Performance.  IEEE Computer, 1994.
  * [29] A. Ketan Shah and M. D. Matani.  An O(1) Algorithm for Implementing the LFU Cache Eviction Scheme.  Technical report, 2010.  ”http://dhruvbird.com/lfu.pdf”.
  * [30] T. Maier, P. Sanders, and R. Dementiev.  Concurrent hash tables: Fast and general(?)!  ACM Trans. Parallel Comput., 5(4):16:1–16:32, 2019.
  * [31] B. Manes.  Caffeine: A High Performance Caching Library for Java 8.  https://github.com/ben-manes/caffeine, 2017.
  * [32] B. Manes.  Segmented Caffeine – Private Communications.  2020\.
  * [33] N. Megiddo and D. S. Modha.  ARC: A Self-Tuning, Low Overhead Replacement Cache.  In FAST, pages 115–130, 2003.
  * [34] A. Metwally, D. Agrawal, and A. E. Abbadi.  Efficient Computation of Frequent and Top-k Elements in Data Streams.  In ICDT, 2005.
  * [35] M. M. Michael.  High performance dynamic lock-free hash tables and list-based sets.  In SPAA, pages 73–82, 2002.
  * [36] M. Mitzenmacher and E. Upfal.  Probability and Computing: Randomized Algorithms and Probabilistic Analysis.  Cambridge University Press, New York, NY, USA, 2005.
  * [37] E. J. O’Neil, P. E. O’Neil, and G. Weikum.  The LRU-K Page Replacement Algorithm for Database Disk Buffering.  ACM SIGMOD Rec., 22(2):297–306, June 1993.
  * [38] S. Park and C. Park.  FRD: A Filtering based Buffer Cache Algorithm that Considers both Frequency and Reuse Distance.  In MSST, 2017.
  * [39] S. Podlipnig and L. Böszörmenyi.  A survey of web cache replacement strategies.  ACM Comput. Surv., 35(4):374–398, Dec. 2003.
  * [40] M. K. Qureshi, D. Thompson, and Y. N. Patt.  The V-Way Cache: Demand-Based Associativity via Global Replacement.  In ISCA, pages 544–555, June 2005.
  * [41] G. Tewari and K. Hazelwood.  Adaptive web proxy caching algorithms.  Technical Report TR-13-04, ”Harvard University”, 2004.
  * [42] V. Tolstopyatov.  xxHash.  https://github.com/OpenHFT/Zero-Allocation-Hashing.
  * [43] G. Urdaneta, G. Pierre, and M. Van Steen.  Wikipedia Workload Analysis for Decentralized Hosting.  Computer Networks, 53(11):1830–1845, 2009.
  * [44] G. Weikum.  UMassTrace Repository.  http://traces. cs. umass. edu/index. php/Storage, 2011.


# Primordial Obliquities of Brown Dwarfs and Super-Jupiters from Fragmenting
Gravito-Turbulent Discs

R. Michael Jennings1, Eugene Chiang1,2
1Department of Astronomy, University of California, Berkeley, Berkeley, CA
94720
2Department of Earth and Planetary Science, University of California,
Berkeley, Berkeley, CA 94720
E-mail<EMAIL_ADDRESS>0000-0002-3959-6572 0000-0002-6246-2310

###### Abstract

Super-Jupiters, brown dwarfs, and stars can form from the collapse of self-
gravitating discs. Such discs are turbulent, with flocculent spiral arms
accelerating gas to transonic speeds horizontally and vertically. Objects that
fragment from gravito-turbulent discs should spin with a wide range of
directions, reflecting the random orientations of their parent eddies. We show
by direct numerical simulation that obliquities of newly collapsed fragments
can range up to 45∘. Subsequent collisions between fragments can further alter
the obliquity distribution, up to 90∘ or down to near-zero. The large
obliquities of newly discovered super-Jupiters on wide orbits around young
stars may be gravito-turbulent in origin. Obliquely spinning fragments are
born on orbits that may be inclined relative to their parent discs by up to
20∘, and gravitationally stir leftover material to many times the pre-
fragmentation disc thickness.

###### keywords:

protoplanetary discs – planets and satellites: formation – turbulence

††pubyear: 2021††pagerange: Primordial Obliquities of Brown Dwarfs and Super-
Jupiters from Fragmenting Gravito-Turbulent Discs–Primordial Obliquities of
Brown Dwarfs and Super-Jupiters from Fragmenting Gravito-Turbulent Discs

##  1 Introduction

Planets that form from unwarped razor-thin circumstellar discs have no choice
but to have zero obliquity at birth, their spin axes and orbit normals
perfectly aligned.111The vector spin angular frequency of a fluid element
equals half its vorticity (e.g. Shapiro 1961). If velocities are purely
horizontal, the vector vorticity is vertical. In Rayleigh-stable discs where
the specific angular momentum increases radially outward, the vorticity points
in the same direction as the orbital angular momentum, i.e. fluid elements
spin prograde. Small obliquities $\ll 1$ radian also result from finite-
thickness discs if the vertical velocities of accreted material sum to near-
zero. The larger the spatial scales over which an object accumulates its mass,
the greater the degree of cancellation.

Planets derived from material on small scales, in discs with significant
vertical motions, may form with large obliquities. A 3D turbulent disc is
filled with eddies whirling in random directions. A tilted eddy that
gravitationally collapses on itself may retain its vector spin angular
momentum to become an obliquely rotating protoplanet. Density fluctuations in
the disc may also torque spins.

Turbulence and gravitational fragmentation go hand in hand. The faster self-
gravitating discs cool, the more they are roiled by spiral density waves. A
gravito-turbulent disc in 3D exhibits transonic radial and vertical motions
extending over many vertical scale heights (Shi & Chiang, 2014). If the disc’s
cooling time becomes shorter than about an orbit time, the disc fragments
(Gammie, 2001). On larger molecular cloud scales, turbulence leads to random
fluctuations in the vector angular momentum of material that accretes to form
a star and its attendant disc, potentially yielding a non-zero stellar
obliquity (e.g. Bate et al., 2010; Thies et al., 2011; Fielding et al., 2015).

Self-gravitating discs are thought to spawn brown dwarf and super-Jupiter
companions to stars on wide orbits. The top-heavy mass function of brown dwarf
secondaries; their tendency to be found at large distances from their host
stars; and their lack of preference for host stellar type all point to top-
down formation by disc fragmentation (e.g. Nielsen et al. 2019; see also
Kratter & Lodato 2016). Recent observations of 2M0122b—a 12–27 Jupiter mass
object orbiting an M dwarf at a projected separation of 50 au—indicate its
obliquity may be large, although the data for this system are not conclusive
because of sky-projection degeneracies (Bryan et al., 2020). The case of HD
106906 is more secure; the F star binary hosts, at separations of hundreds of
au’s, a $\sim$12 Jupiter mass companion viewed nearly pole-on, and a debris
disc viewed nearly edge-on (Bryan et al. 2021, in press). Large obliquities of
super-Jupiters on wide orbits suggest fragmentation of a gravito-turbulent
disc.

We aim here to provide further evidence supporting this interpretation. We
present 3D numerical fluid simulations of secularly cooling, gravito-turbulent
discs and measure the obliquities of collapsed fragments. A shearing box is
employed to resolve the turbulent structure of a local patch of the disc. The
drawback of the shearing box is that it cannot track the evolution of fragment
orbits; orbital radius and by extension orbital angular momentum are not
defined in a shearing box formalism which neglects radial curvature terms.
Thus for obliquity, which measures the angle between the spin and orbital
angular momentum vectors, we can only track the spin vector; the orbit normal
is assumed fixed along the vertical axis of the box. This shortcoming is
removed in a global calculation at the cost of spatial resolution. Hall et al.
(2017) presented global 3D smooth-particle-hydrodynamics simulations of self-
gravitating discs, apparently finding all bound fragments to form with
obliquities $\lesssim 10^{\circ}$. Collisions, and mass and angular momentum
transfer between fragments led to four fragments among the forty recorded
spinning retrograde with obliquities $>160^{\circ}$. One fragment was plotted
to have an obliquity of $35^{\circ}$, but its dynamical history was not
discussed. We offer here a complementary calculation that more finely resolves
the 3D velocity fields in self-gravitating discs, in particular the vertical
flows that can lead to oblique rotators.

The rest of this paper is organized as follows. Section 2 describes our
numerical code and setup, section 3 presents our results on obliquities and
other fragment properties, and section 4 provides a summary and outlook.

##  2 Code and Setup

Our computational tools and procedures are similar to those of Shi & Chiang
(2014), and we provide here just a brief description highlighting differences.
We use the code Athena to solve the hydrodynamic equations in a 3D shearing
box, including vertical tidal gravity from the central star (Stone et al.,
2008). We adopt the van Leer integrator (van Leer, 2003; Stone & Gardiner,
2009); a piecewise linear spatial reconstruction in the primitive variables;
the Harten–Lax–van Leer contact (HLLC) Riemann solver; and orbital advection
algorithms that speed runtimes and improve conservation (Masset, 2000; Johnson
et al., 2008; Stone & Gardiner, 2010). Disc self-gravity is included by
solving Poisson’s equation using fast Fourier transforms (Koyama & Ostriker,
2009; Kim et al., 2011). Boundary conditions are shearing-periodic in radial
displacement $x$, and periodic in azimuth $y$. In vertical height $z$, vacuum
boundary conditions are adopted for the self-gravitational potential $\Phi$,
and reflecting boundary conditions for the other hydrodynamic flow variables.
This last choice differs from Shi & Chiang (2014) who used vertical outflow
boundary conditions. We found that outflow conditions resulted in significant
mass loss from the simulation domain after fragments form and gravitationally
stir the remaining disc material (see section 3.4). Mass loss was not an issue
for Shi & Chiang (2014) whose focus was on turbulence in non-fragmenting
discs. With reflecting vertical boundary conditions for density and velocity,
the mass in the box is conserved to a few percent.

Initial conditions at time $t=0$ are of a horizontally uniform disc with a
vertical density profile in hydrostatic equilibrium including self-gravity:

$\frac{1}{\rho}\frac{dP}{dz}=-\Omega^{2}z-4\pi
G\int^{z}_{0}\rho(z^{\prime})dz^{\prime}$ (1)

for density $\rho$, pressure $P$, fixed orbital angular frequency $\Omega$,
and gravitational constant $G$. Given an initial surface density $\Sigma_{0}$
and initial midplane density $\rho_{0}$, we define a characteristic initial
disc half-thickness $H\equiv\Sigma_{0}/(2\rho_{0})$ over which most of the
disc mass is contained. The gas is assumed to obey a polytrope
$P=K\rho^{\gamma}$ initially, with an initial midplane sound speed $c_{\rm
s0}=\sqrt{\gamma K\rho_{0}^{\gamma-1}}$. Upon setting $\gamma=5/3$ and
$Q_{0}\equiv c_{\rm s0}\Omega/(\pi G\Sigma_{0})=1$, equation (1) may be solved
for the initial density profile $\rho(z)$ (Shi & Chiang, 2013, 2014). The
solution is plotted as a dotted curve in Figure 1. Analytically, $\rho$
vanishes at a finite height $z\simeq\pm 2H$. To avoid having to resolve steep
density and pressure gradients near the vertical boundaries of the box, we
introduce a density floor of $\rho_{\rm min}=10^{-4}\rho_{0}$. We work in code
units for which $\rho_{0}=\Omega=H=1$; it follows that $\Sigma_{0}=2$, $c_{\rm
s0}=2.126$, $K=2.712$, and $G=0.3384$. To seed gravitational instability we
introduce (at $t=0$ only) random cell-to-cell velocity perturbations with a
maximum amplitude of either $0.02$ or $0.2$$c_{\rm s0}$ for $|z|<2H$; we have
verified in a few cases that our results do not depend on the exact choice of
initial perturbation.

Our fiducial box size is $\\{L_{x},L_{y},L_{z}\\}=\\{16H,64H,24H\\}$. We want
our box to be large enough to avoid much disc material, and most especially
bound fragments, crashing into the reflecting box lids, and to accommodate as
many turbulent sub-structures (potential fragments) as possible. If we say
that our simulation is intended to model a portion of the disc at radius $a$
having a typical protoplanetary disc aspect ratio $H/a=0.1$, then our box
extends the full disc circumference in azimuth, and $0.2a$–$1.8a$ in radius
(in order-unity violation of the shearing sheet formalism which assumes the
computational domain is less than the disc radius in size). The azimuthal
length $L_{y}$ is large enough to accommodate the long-wavelength modes
responsible for maintaining steady gravito-turbulence without “bursty”
behaviour (Booth & Clarke, 2019). Our fiducial resolution is 4 cells per $H$;
we also experiment with high-resolution runs having 8 cells per $H$. Shi &
Chiang (2014) used identical resolutions and found no appreciable difference
between them; similar results were obtained by Booth & Clarke (2019) who
tested up to 32 cells per $H$.

In our fiducial runs, cooling is effected by lowering the internal energy
density $U=P/(\gamma-1)$ of every cell at a rate $\dot{U}=-U/t_{\rm cool}$ for
cooling time $t_{\rm cool}\equiv\beta/\Omega$. In addition to the
aforementioned floor on mass density $\rho$, we impose a floor on the energy
density (and by extension pressure) of $U_{\rm min}=10^{-4}$ in code units. We
start runs with $\beta=\beta_{0}\in\\{10,20\\}$ for times $t<40\Omega^{-1}$;
during this earliest phase of the simulation, the disc is in a quasi-steady,
non-fragmenting, gravito-turbulent state. To induce fragmentation, we then
accelerate cooling by lowering $\beta$ from $\beta_{0}$ to a minimum value of
2 over a linear ramp in time from $t=40\Omega^{-1}$ to $t=120\Omega^{-1}$.
Typically one or two fragments form when $\beta=2$.

We also experiment with an alternative cooling prescription introduced by Shi
& Chiang (2014) to simulate radiative cooling in an optically thin medium with
constant opacity:

$\dot{U}=-\frac{1}{b(\gamma-1)}\frac{P^{4}}{\rho^{3}}\,.$ (2)

The cooling time $U/|\dot{U}|$ then scales as
$[P/(\gamma-1)]/|\dot{U}|=b(\rho/P)^{3}\propto b/T^{3}$ for coefficient $b$
and temperature $T$; this scaling is appropriate for gas whose energy density
scales as $\rho T$ and whose volume emissivity scales as $\rho T^{4}$. We
start with $b=b_{0}\in\\{500,600\\}$ for the first $40\Omega^{-1}$, and then
accelerate cooling by ramping $b$ down at a rate of
$db/dt=100/(80\Omega^{-1})$. Note that these optically thin cooling
experiments are performed to explore an issue of principle, namely the
sensitivity of obliquity outcomes to the cooling prescription, and not because
we think massive, self-gravitating discs are necessarily optically thin in
practice. Early in the evolution of a disc, when it is still being fed by a
natal envelope, we suspect it is optically thick to its own cooling radiation.

Figure 1: Vertical density profiles of the background circumstellar (CS) disc
in our high-resolution (8 cells per $H$), $\beta_{0}=20$ run. The initial
condition (vertical hydrostatic equilibrium) is plotted as a dotted line. When
the disc becomes gravito-turbulent, the disc thickens. Thin colored lines are
horizontally averaged snapshots plotted before the disc fragments; as
$\beta=\Omega t_{\rm cool}$ decreases from 20 to 2 and cooling accelerates,
the disc flattens. The $\beta=2$ curve describes the vertical density profile
just before the fragment forms. At this stage the disc is relatively thin and
the density profile within $z=\pm 2H$ resembles the hydrostatic initial
condition. Once a fragment forms the disc puffs back up, and its midplane
oscillates about $z=0$, as seen by the dashed lines representing horizontally
averaged snapshots with the fragment removed. The stirring is caused by the
gravity of the fragment.

##  3 RESULTS

Figure 2: Column density snapshots of our high resolution (8 cells per $H$),
$\beta_{0}=20$ run. The top row, with snapshot times listed in the header and
corresponding $\beta=\Omega t_{\rm cool}$ values annotated in the corners,
shows edge-on views of the disc in the meridional $x$-$z$ plane, with density
integrated in the azimuthal-$y$ direction. The bottom row shows face-on views
of the disc, plotting surface density in $x$-$y$ plane (density integrated
along the vertical $z$ direction) for subsections of the simulation domain. At
early times, trailing spiral arms amplify. Later, a single fragment forms
whose spin angular momentum vector $\bm{J}$ points $\sim$35 degrees from the
vertical. The white arrow of arbitrary length shows the $x$-$z$ projection of
$\bm{J}$ and the dashed circle encloses the fragment’s Hill sphere. The
fragment is disc-shaped (top row), exhibits two tidal arms (bottom row), and
has a spin axis which is directed in the positive $z$-direction (the spin is
prograde; see Fig. 3) and which precesses about the $z$-axis with angular
frequency $\Omega$ in the retrograde direction. The precession causes the
apparent orientation of the fragment (a.k.a. circumplanetary disc) to change
between the third and fourth columns and is a purely kinematic effect arising
from the rotating frame. Figure 3: Left: Iso-contours of the self-
gravitational potential $\Phi$, taken in an $x$-$z$ slice (at fixed $y$
passing through the fragment center of mass; top panel) and an $x$-$y$ slice
(at fixed $z$ passing through the fragment center of mass; bottom panel), for
the high-resolution, $\beta_{0}=10$ run at $\Omega t=112$. Color background
traces column densities. Center: Volume density iso-contours for these same
slices. Color background traces volume density. Right: Velocities for our
slices. The circumplanetary disc (CPD) rotates prograde (bottom panel) and
emits a low-density bipolar outflow (top panel). The bipolar wind appears
driven by steep pressure gradients near CPD surface layers.

###  3.1 Overview and Diagnostics

Simulations typically play out as follows. They begin with a period of
gravito-turbulence which intensifies as $\beta=\Omega t_{\rm cool}$ (or the
analogous parameter $b$ in the optically thin runs) decreases and cooling
accelerates. At some point during the cooling the disc fragments, usually when
$\beta$ reaches its prescribed minimum of 2. The vertical disc velocities
generated during the pre-fragmentation, gravito-turbulent phase will lead to
non-zero fragment obliquities. Shi & Chiang (2014) found that such velocities
approach but remain below the sound speed for most of the disc mass (see their
figure 7). We confirm that the disc just prior to fragmentation is still
relatively vertically thin — see the progression of thin colored lines in Fig.
1 from high to low $\beta$. Snapshots of the evolution of our high-resolution,
$\beta_{0}=20$ run are shown in Fig. 2.

Depending on the run, between one and four fragments initially form. Multiple
fragments subsequently merge. Occasionally a fragment “spins off” another when
one of its extended tidal tails collapses under its own self-gravity.
Fragments are formally located by searching for local minima in the self-
gravitational potential $\Phi$. The outer boundary of a fragment is the
equipotential surface whose largest distance $R_{\rm max}$ from the enclosed
local potential minimum coincides with the Hill radius $R_{\rm Hill}=(GM_{\rm
frag}/3\Omega^{2})^{1/3}$, where $M_{\rm frag}$ is the mass contained within
this outermost equipotential (Fig. 3, left panel). The radii $R_{\rm max}$ and
$R_{\rm Hill}$ are computed iteratively until $R_{\rm max}=R_{\rm Hill}$ to
within 1 simulation cell width. We confirm that the Bondi radius
$R_{\text{Bondi}}=2GM_{\rm frag}/\langle c_{\rm s}^{2}\rangle_{\rho}$, where
$\langle...\rangle_{\rho}$ denotes a density-weighted average over the
fragment, is always larger than $R_{\rm Hill}$; fragment sizes are limited by
stellar tides and not by thermal pressure, as appears clear from the two-armed
tidal tails extending out to the fragment boundary (Figs. 2 and 3). The
fragment center of mass is computed over all cells within the fragment
boundary, which we reiterate is not assumed spherical but conforms to a level
surface of $\Phi$. In practice, the center of mass differs negligibly from the
potential minimum.

The inertial-frame spin angular momentum vector of a fragment is computed as

$\displaystyle\bm{J}=\sum_{i}m_{i}\bm{R}_{i}\times\left(\bm{V}_{i}+\Omega\bm{\hat{z}}\times\bm{R}_{i}\right)$
(3)

where $m_{i}$, $\bm{R}_{i}$, and $\bm{V_{i}}$ are the mass, displacement
vector, and velocity vector of the $i$th cell, measured in the shearing-box
frame relative to the center of mass, and the sum is over all cells within the
fragment boundary. The term $+\Omega\bm{\hat{z}}\times\bm{R}_{i}$, where
$\Omega=1$ is the frame rotation frequency, corrects for the apparent
retrograde “precession” of the fragment’s spin axis seen in the shearing-box
frame (e.g. Nesvorny et al. 2019, but note that their equation 8 has a sign
error for the correction term). This precession, which occurs with angular
frequency $-\Omega\bm{\hat{z}}$, is clearly evident in movies of our
simulations (see Fig. 2 for two snapshots at two different precession phases).
It is a purely kinematic effect that arises because the shearing-box frame
rotates with angular frequency $+\Omega\bm{\hat{z}}$. In practice, removing
the frame rotation to recover the inertial-frame $\bm{J}$ is a small
correction.

The obliquity $\Psi$ of a fragment is the angle between the spin $\bm{J}$ and
the orbit normal. As noted in the Introduction, a local shearing box cannot
track the evolution of the orbital angular momentum vector (there is no unique
orbital radius in a local calculation), and so we have no choice but to assume
the orbit normal always points in the $\bm{\hat{z}}$ direction, as it does by
construction at the beginning of our simulations. Thus we compute the
obliquity as

$\cos{\Psi}=\frac{\bm{J}\cdot\bm{\hat{z}}}{|\bm{J}|}\,.$ (4)

###  3.2 Obliquities

Figure 4: Time evolution of fragment properties. From top to bottom: Hill
radius $R_{\rm Hill}$, fragment mass $M_{\rm frag}$ relative to the total mass
in the simulation box, obliquity $\Psi$, height $z_{\rm frag}$ of the fragment
center of mass, and total spin angular momentum $|\bm{J}|$ normalized to a
rough measure of the spin angular momentum needed for break-up. Different
colors refer to different simulations as annotated. Different linestyles
(solid vs. dashed) refer to different fragments within the same simulation.
For our fiducial run (black), two fragments appear at $\Omega t\sim 125$ and
eventually merge into one at $\Omega t\sim 170$. The optically thin,
$b_{0}=500$ run (yellow) is especially dynamic: two fragments first appear at
$\Omega t\sim 60$ and merge soon thereafter, and the resultant object
periodically spins off and re-accretes a lower-mass companion. Large
obliquities can arise either at the time a fragment forms and accretes (black,
purple, red) or when fragments collide (yellow). The disc in the high-
resolution, $\beta_{0}=10$ run (blue) is initially fast-cooling and fragments
before turbulence has had a chance to develop; the resultant single fragment
has low obliquity. Spikes in the bottommost panel are due either to close
encounters between fragments which rapidly alter spin angular momenta
(yellow), or to sudden decreases and increases in fragment mass which change
the denominator $\propto M_{\rm frag}^{5/3}$ (blue). See text for more
details.

Fig. 4 plots fragment obliquities $\Psi$ from our various simulations. In our
fiducial $\beta$-cooling run (black curve: $L_{z}=24H$, resolution = 4 cells
per $H$, initial $\beta_{0}=10$), the disc initially spawns two fragments
having small obliquities. As the fragments interact with each other and
leftover circumstellar (CS) gas, the obliquities increase from a combination
of gravitational forcing and direct accretion of mass and angular momentum.
Eventually the fragments merge into a single object having $\Psi\simeq
25$–45∘.

In other $\beta$-cooling runs, the disc fragments into a single object whose
obliquity appears determined from the spin angular momentum of accreted gas
and not from gravitational forcing by ambient material. Obliquities on the
order of 1 radian are produced in both the small box run (red curve:
$L_{z}=12H$), and the high-resolution run with slower initial cooling (dashed
blue: $8$ cells per $H$, $\beta_{0}=20$).

Fragment masses $M_{\rm frag}$ toward the end of the aforementioned
simulations range from 30% to 65% of the total (conserved) box mass $M_{\rm
total}$, fluctuating as fragments repeatedly collide with disc filaments.
Sometimes the colliding filament arises from the fragment itself: because of
our shearing-periodic radial boundaries, a fragment’s tidal arm can extend
past one radial boundary to re-emerge from the other boundary and swipe the
original fragment. This self-swiping is not necessarily unphysical; the arm
that collides with the fragment could be interpreted as originating from
another fragment outside the simulation domain. In any case, the fluctuations
in mass are confined to the periphery of the fragment and do not appear to
much affect the obliquity (see also the next section 3.3 on fragment internal
structure).

In the high-resolution $\beta$-cooling run with faster initial cooling (solid
blue: $8$ cells per $H$, $\beta_{0}=10$), $\Psi$ is small, $\lesssim
10^{\circ}$, a consequence of the fragment forming early, before gravito-
turbulent vertical motions fully develop (see also Fig. 5). This case
underscores how the approach to fragmentation matters: large obliquities more
easily arise in discs that are already turbulent before they collapse. The
obliquity remains small despite the fragment losing $\sim$90% of its mass in
an apparent self-swipe; the fragment’s small obliquity and close confinement
to the disc midplane (see the panel in Fig. 4 plotting $z_{\rm frag}$, the
height of the fragment center of mass) make it especially susceptible to self-
swiping. The dynamics here are essentially of 2D midplane flows. During the
mass loss and subsequent re-accretion, the obliquity varies between zero and
$10^{\circ}$.

Discs in the optically thin $b$-experiments fragment relatively easily, before
vertical turbulent motions develop. Hot, overpressured gas created at the
intersection of converging flows and colliding filaments does not stay hot for
long, as hotter gas cools proportionately faster according to our cooling
prescription (2). In a sense, the equation of state is soft, and so when disc
gas compresses in one direction, it is more prone to gravitationally collapse
than to flow in another direction. Initially in these optically thin runs,
midplane gas collapses into more numerous and lower-mass fragments, all having
small obliquities, as compared to our $\beta$-cooling models. In the
$b_{0}=600$ run, between one and four low-obliquity fragments are present near
the midplane at any given time, with the number changing from mergers and
tidal spin-offs (Fig. 5, middle column). In the $b_{0}=500$ run, two fragments
initially form with low obliquities, and soon collide to form an object with a
near-90∘ obliquity. The subsequent evolution is chaotic (Figs. 4 and 5), with
multiple spin-offs and mergers leading eventually to a single low-obliquity
object.

###  3.3 Fragment Internal Structure

Fig. 6 profiles the radial structure of a fragment (the same one displayed in
Fig. 2), plotting the mass, scalar spin angular momentum, and obliquity of
material enclosed within equipotentials of varying size $R$ (the largest
distance between a given equipotential surface and the center of mass).
Roughly half of the total fragment mass is contained within an “inner core” of
size $R\simeq 0.2R_{\rm Hill}$. This inner core contains 20–25% of the total
spin angular momentum, and its mass profile $M_{\rm frag}(R)$ appears to scale
linearly with $R$, suggesting the core is supported against self-gravity in
large part by thermal pressure (an isothermal sphere in hydrostatic
equilibrium has an enclosed mass that scales linearly with radius).

The other half of the fragment mass is located in the “outer circumplanetary
disc (CPD)” which extends from $R\simeq 0.2R_{\rm Hill}$ to $R_{\rm
max}=R_{\rm Hill}$. About 80% of the total spin angular momentum (whose
magnitude $|\bm{J}|$ is comparable to the break-up value; see Fig. 4) is
carried by the outer CPD. There is a slight warp across the outer CPD as
$\Psi(R)$ decreases by several degrees from small to large $R$; see also the
volume rendering in Fig. 7.

###  3.4 Vertical Dynamics

Disc vertical velocities arising from gravito-turbulence result in fragments
that not only rotate obliquely but also have vertical center-of-mass motions.
From Fig. 4, in those runs where $\Psi\sim 1$ rad, the vertical position of a
fragment $z_{\rm frag}$ oscillates with a typical amplitude of $\sim$1–$2H$;
in one high-resolution run, vertical excursions extend up to $\sim$$4H$. In
the runs where $\Psi<10^{\circ}$, vertical oscillations are correspondingly
small, $|z_{\rm frag}|\ll H$.

Vertical oscillations imply that fragments have non-zero orbital
inclinations—unsurprisingly, given their formation within a turbulent disc.
The orbital inclinations cannot be read directly from our local shearing-box
simulations which have no well-defined orbital radius $a$. However, if we
assume an aspect ratio of $H/a\sim 0.1$ typical of protoplanetary discs, then
vertical oscillation amplitudes of 1–4$H$ imply orbital inclinations of
0.1–0.4 rad $\sim$ 6–20∘. In making this estimate we have assumed that the
midplane pressure scale height $|dz/d\ln P|\sim 0.1a$ of an actual disc can be
approximated by our code-based $H\equiv\Sigma_{0}/(2\rho_{0})$, where the
latter is defined as part of our hydrostatic initial condition. This
identification is good if real discs have density profiles resembling our
assumed initial condition within a few pressure scale heights of the midplane
where most of the disc mass resides. Fig. 1 lends support to this assumption
by showing that fast-cooling, gravito-turbulent, pre-fragmentation discs have
vertical density profiles within $z=\pm 2H$ resembling those of our
hydrostatic initial condition.

As obliquity is the angular separation between an object’s spin axis and orbit
normal, it obviously depends on orbital inclination. A shortcoming of our
formula (4) for $\Psi$ imposed by the shearing box that it does not account
for non-zero inclination. We expect, however, the error to be small to the
extent that our estimated inclinations of 6–20∘ are smaller than our measured
45–$90^{\circ}$ spin axis tilts. Furthermore, our estimated inclinations
should not bias our measured values for $\Psi$ upward or downward, as the
phase of the vertical oscillations in $z_{\rm frag}$ does not correlate with
the “precession” phase of the spin axis (see section 3.1).

Fragments gravitationally stir leftover circumstellar (CS) disc material
strongly. For the high-resolution, $\beta_{0}=20$ simulation plotted in Fig.
1, after the fragment forms, the ambient CS disc thickens dramatically; gas
densities at the vertical boundaries of our simulation domain remain within a
factor of $\sim$10 of the midplane density. The CS disc midplane itself is
forced out of the $z=0$ plane; Fig. 1 shows the maxima of the post-
fragmentation vertical density profiles (dashed lines) varying between
$z\simeq\pm 2H$. The gravity of the fragment is responsible for this forcing
of CS disc gas. The fragment mass is an order-unity fraction of the original
CS disc mass (Fig. 4), and its surface escape velocity $\sim$$\sqrt{2GM_{\rm
frag}/R_{\rm Hill}}$ is $\sim$$2.5\times$ larger than the CS disc sound speed
evaluated at the midplane just prior to fragmentation.

We emphasize that the most extreme thickening of the disc, when an order-unity
fraction of the disc mass is thrown to heights exceeding several $H$, occurs
only after the fragment forms; this violent puffing-up is a consequence and
not a cause of fragmentation, and does not affect the fragment obliquity. For
example, the high-resolution $\beta_{0}=20$ run which yields high obliquity
(Fig. 2) and the high-resolution $\beta_{0}=10$ run which produces low
obliquity (Fig. 5) exhibit post-fragmentation disc thickening of similar
magnitude. Obliquities derive instead from whatever vertical velocities are
present in the disc pre-fragmentation. Near the midplane, within $z=\pm 2H$,
these velocities approach but remain less than the sound speed (Shi & Chiang
2014, their figure 7). Fig. 1 shows that the disc just prior to fragmentation
is much thinner than the disc after fragmentation — contrast the thin yellow
lines for $\beta\approx 2$ with the thick dashed curves. Just before it
fragments, the disc between $z=\pm 2H$ still contains the lion’s share of the
mass, and its vertical density profile is similar to that of the hydrostatic
initial condition. Non-zero obliquity follows less from puffing up — there
isn’t much puffing before fragmentation — and more from the disc acquiring
non-zero turbulent vertical velocities because of self-gravity.

Figure 5: More column density snapshots like those in Fig. 2, using the same
color scale. Obliquities are small in the high-resolution $\beta_{0}=10$ run
and in the optically thin $b_{0}=600$ run (for the latter we show only 3 of
the 4 fragments that form initially). In the optically thin $b_{0}=500$ run,
several fragments momentarily coalesce into a rotating ring having an
obliquity of nearly 90∘. Figure 6: Fragment internal structure, plotting the
mass, spin angular momentum, and obliquity of material enclosed with an
equipotential of running size $R$. More than half of the total fragment mass
and about half of the total spin angular momentum are contained within $R\sim
0.3R_{\rm Hill}$. Data are shown for the high-resolution, $\beta_{0}=20$ run
at two times as indicated. There are about 32 simulation cells spanning
$R_{\rm Hill}\simeq 4H$. Figure 7: Volume rendering of the fragment in our
high resolution, $\beta_{0}=20$ run at $\Omega t=191$. Gas with
$$1.8\text{\times}{10}^{3}$\geq\rho/\rho_{0}\geq 10^{2}$ is colored yellow,
$10^{2}>\rho/\rho_{0}\geq 1$ is green, and $1>\rho/\rho_{0}\geq 10^{-2}$ is
blue. The yellow inner core is $\sim$$1H$ in diameter.

##  4 Summary and Discussion

Brown dwarf and super-Jupiter companions to stars are thought to form by
gravitational instability in circumstellar discs (e.g. Nielsen et al. 2019).
Self-gravitating discs are turbulent, containing a “swirling hotch-potch of
pieces of spiral arms” (Goldreich & Lynden-Bell, 1965; Gammie, 2001; Shi &
Chiang, 2014). We have shown by direct numerical simulation that gravito-
turbulence can lead to gravitationally collapsed fragments that spin with
order-unity obliquities. Initial obliquities are determined by the spin
angular momentum of accreted material; a nascent fragment spins in the same
direction as the turbulent whorl from which it was born. In our simulations,
this primordial spin can point as much as $45^{\circ}$ from the orbit normal.
Subsequent close encounters, collisions, and mergers of fragments can increase
or decrease obliquities, up to $\sim$$90^{\circ}$ or down to near-zero.

The caveat to these results is that they were obtained in a shearing box for
which stellocentric radius and by extension orbital angular momentum and
orbital inclination are not defined. Thus we have assumed in computing
obliquities that orbit normals are fixed in the direction they were pointing
at the start of our simulations, parallel to the box’s vertical dimension. At
the same time we have seen both the fragments and the leftover circumstellar
disc gas oscillate vertically in our simulations, implying that non-zero
orbital inclinations are excited in the course of disc evolution. These
vertical motions are expected from gravito-turbulence and gravitational
forcing by fragments. We have measured the vertical oscillations to have an
amplitude comparable to the vertical disc scale height $H$. Insofar as real-
life protoplanetary discs have height-to-radius ratios of $H/a\sim 0.1$, we
estimate that the corresponding orbital inclinations are on the order of 0.1
rad, somewhat smaller than the order-1 spin tilts measured in several of our
runs. Furthermore, these vertical oscillations appear randomly phased with
respect to spin axis orientations, implying that our turbulence-induced
obliquities are not systematically reduced by non-zero orbital inclinations.
Thus our conclusion that order-unity obliquities can result from gravitational
instability seems secure, pending global 3D simulations that can
simultaneously track spin and orbit. Note further that in all our runs,
fragments move radially, sometimes across the entire simulation domain
spanning $16H$; these radial motions are only to be expected from the
background turbulence, and suggest that objects formed by gravitational
instability should have large orbital eccentricities.

The approach to fragmentation, i.e. the exact manner in which Toomre’s $Q$ and
the gas cooling time $t_{\rm cool}$ drop to the critical values needed for
gravitational collapse, matters for obliquity. If the descent to instability
is fast (lasting less than ten or so orbital periods), or if overdensities
cool so efficiently that gas pressure cannot compete against self-gravity (in
other words if the equation of state is too soft), then vertical turbulent
motions may not develop in the disc before gas collapses, more-or-less
symmetrically, into a low-obliquity fragment. The evolution to fragmentation
depends on the radiation-thermodynamics of how gas cools, how mass is
transported across disc radius, and ongoing accretion from the disc’s natal
envelope. These issues are beyond the scope of our study. Kratter et al.
(2010, their section 6) discuss a fragmentation scenario whereby a disc is fed
at too high a rate by an external infalling envelope for disc spiral density
waves to torque the mass away. The local surface density then builds until $Q$
falls to order-unity; meanwhile, the cooling criterion is assumed to be
readily satisfied (their figure 1; see also Tobin et al. 2016, their extended
data figure 4). The disc becomes gravito-turbulent, and potentially spawns an
obliquely rotating fragment, if the timescale for the local disk to grown in
mass and $Q$ to approach unity is at least as long as the local orbital time;
this seems a realistic prospect. Note that this scenario, which assumes the
cooling criterion is always satisfied while $Q$ decreases, is the converse of
what we implemented in our code for numerical convenience; the underlying
point is that disc evolutionary timescales can be slow enough to trigger
turbulence before fragmentation.

Gravito-turbulence sets an object’s initial obliquity. A litany of other
effects may subsequently alter its value, either by changing the spin axis or
the orbital axis. These include, e.g., torques from the host star that can
either damp or excite the inclination of a circumplanetary disc relative to
the circumstellar orbit plane (Lubow & Ogilvie, 2000; Martin & Armitage,
2021); collisions between protoplanets (e.g. Li et al. 2020; Li & Lai 2020);
changes to the orbit plane from ongoing circumstellar disc accretion
(Tremaine, 1991); and spin-orbit resonances brought about by migration
(Millholland & Laughlin, 2019) or planet-planet scattering (Hong et al., 2021;
Li, 2021). Planetary magnetic fields are thought to be critical in regulating
the spin rates of proto-gas-giants and thereby enabling their contraction to
observed sizes (Ginzburg & Chiang, 2020); magnetic torques may change the spin
direction as well (e.g. Spalding & Batygin 2015).

Obliquities of objects having masses upwards of $\sim$10 Jupiter masses are
beginning to be constrained observationally. The first brown dwarfs/super-
Jupiters to be examined for their spin-orbit angles are 2M0122b (Bryan et al.,
2020) and HD 106906b (Bryan et al. 2021, in press). Both companions orbit
their host stars at the large, $>50$ AU distances conducive to disc
fragmentation, and both evince, modulo large uncertainties, order-unity
obliquities. In the case of the HD 106906 system, it has been established with
high confidence that the companion is viewed nearly pole-on (Bryan et al.
2021, in press) while a debris disc orbiting the same host is viewed nearly
edge-on (Kalas et al., 2015; Lagrange et al., 2016; Nguyen et al., 2021). Disc
gravitational instability is also relevant for the formation of binary stars.
Whether stellar-mass secondaries should have large obliquities is unclear; it
may be that such companions accrete so much of the parent disc that turbulent
velocity fluctuations in the accreted material average to zero. Observational
errors on the degree of spin-orbit alignment in solar-type binaries are
currently too large to draw conclusions (Justesen & Albrecht, 2020).

## Acknowledgements

We thank Ji-Ming Shi for providing a copy of his code, and Marta Bryan, Dan
Fabrycky, Sivan Ginzburg, Greg Laughlin, Dong Lai, Gongjie Li, Chris Matzner,
and Zhaohuan Zhu for helpful discussions. Kaitlin Kratter provided a
thoughtful referee’s report which improved the content of this paper. The
simulations were run on the Savio computer cluster provided by the Berkeley
Research Computing program, which is supported by the UC Berkeley Chancellor,
Vice Chancellor for Research, and Chief Information Officer. This research
made use of the open source project yt (yt-project.org; Turk et al. 2011), and
was supported by NASA grant NNX15AD95G/NEXSS.

## Data Availability

Movies of the simulations may be accessed at
https://michaeljennings11.github.io/pages/gi.

## References

  * Bate et al. (2010) Bate M. R., Lodato G., Pringle J. E., 2010, MNRAS, 401, 1505
  * Booth & Clarke (2019) Booth R. A., Clarke C. J., 2019, Monthly Notices of the Royal Astronomical Society, 483, 3718
  * Bryan et al. (2020) Bryan M. L., et al., 2020, The Astronomical Journal, 159, 181
  * Fielding et al. (2015) Fielding D. B., McKee C. F., Socrates A., Cunningham A. J., Klein R. I., 2015, MNRAS, 450, 3306
  * Gammie (2001) Gammie C. F., 2001, The Astrophysical Journal, 553, 174
  * Ginzburg & Chiang (2020) Ginzburg S., Chiang E., 2020, Monthly Notices of the Royal Astronomical Society: Letters, 491, L34
  * Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, Monthly Notices of the Royal Astronomical Society, 130, 125
  * Hall et al. (2017) Hall C., Forgan D., Rice K., 2017, Monthly Notices of the Royal Astronomical Society, 470, 2517
  * Hong et al. (2021) Hong Y.-C., Lai D., Lunine J. I., Nicholson P. D., 2021, ApJ, submitted, arXiv:2103.15902
  * Johnson et al. (2008) Johnson B. M., Guan X., Gammie C. F., 2008, The Astrophysical Journal Supplement Series, 177, 373
  * Justesen & Albrecht (2020) Justesen A. B., Albrecht S., 2020, A&A, 642, A212
  * Kalas et al. (2015) Kalas P. G., et al., 2015, ApJ, 814, 32
  * Kim et al. (2011) Kim C.-G., Kim W.-T., Ostriker E. C., 2011, ApJ, 743, 25
  * Koyama & Ostriker (2009) Koyama H., Ostriker E. C., 2009, Astrophysical Journal, 693, 1316
  * Kratter & Lodato (2016) Kratter K., Lodato G., 2016, Annual Review of Astronomy and Astrophysics, 54, 271
  * Kratter et al. (2010) Kratter K. M., Murray-Clay R. A., Youdin A. N., 2010, ApJ, 710, 1375
  * Lagrange et al. (2016) Lagrange A.-M., et al., 2016, Astronomy & Astrophysics, 586, L8
  * Li (2021) Li G., 2021, ApJ, submitted, arXiv:2103.15843
  * Li & Lai (2020) Li J., Lai D., 2020, ApJ, 898, L20
  * Li et al. (2020) Li J., Lai D., Anderson K. R., Pu B., 2020, Monthly Notices of the Royal Astronomical Society, 501, 1621
  * Lubow & Ogilvie (2000) Lubow S. H., Ogilvie G. I., 2000, The Astrophysical Journal, 538, 326
  * Martin & Armitage (2021) Martin R. G., Armitage P. J., 2021, The Astrophysical Journal Letters, 912, L16
  * Masset (2000) Masset F., 2000, Astronomy and Astrophysics Supplement Series, 141, 165
  * Millholland & Laughlin (2019) Millholland S., Laughlin G., 2019, Nature Astronomy, 3, 424
  * Nguyen et al. (2021) Nguyen M. M., De Rosa R. J., Kalas P., 2021, AJ, 161, 22
  * Nielsen et al. (2019) Nielsen E. L., et al., 2019, AJ, 158, 13
  * Shapiro (1961) Shapiro A., 1961, Illustrated Experiments in Fluid Mechanics. National Committee for Fluid Mechanics Films
  * Shi & Chiang (2013) Shi J.-M., Chiang E., 2013, ApJ, 764, 20
  * Shi & Chiang (2014) Shi J.-M., Chiang E., 2014, ApJ, 789, 34
  * Spalding & Batygin (2015) Spalding C., Batygin K., 2015, Astrophysical Journal, 811, 82
  * Stone & Gardiner (2009) Stone J. M., Gardiner T., 2009, New Astronomy, 14, 139
  * Stone & Gardiner (2010) Stone J. M., Gardiner T. A., 2010, Astrophysical Journal, Supplement Series, 189, 142
  * Stone et al. (2008) Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, The Astrophysical Journal Supplement Series, 178, 137
  * Thies et al. (2011) Thies I., Kroupa P., Goodwin S. P., Stamatellos D., Whitworth A. P., 2011, MNRAS, 417, 1817
  * Tobin et al. (2016) Tobin J. J., et al., 2016, Nature, 538, 483
  * Tremaine (1991) Tremaine S., 1991, Icarus, 89, 85
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  * van Leer (2003) van Leer B., 2003, 16th AIAA Computational Fluid Dynamics Conference, 1, 192