## Abstract

Dissipation is ubiquitous in nature and plays a crucial role in quantum systems such as causing decoherence of quantum states. Recently, much attention has been paid to an intriguing possibility of dissipation as an efficient tool for the preparation and manipulation of quantum states. We report the realization of successful demonstration of a novel role of dissipation in a quantum phase transition using cold atoms. We realize an engineered dissipative Bose-Hubbard system by introducing a controllable strength of two-body inelastic collision via photoassociation for ultracold bosons in a three-dimensional optical lattice. In the dynamics subjected to a slow ramp-down of the optical lattice, we find that strong on-site dissipation favors the Mott insulating state: The melting of the Mott insulator is delayed, and the growth of the phase coherence is suppressed. The controllability of the dissipation is highlighted by quenching the dissipation, providing a novel method for investigating a quantum many-body state and its nonequilibrium dynamics.

## INTRODUCTION

Dissipation—coupling to the environment—plays an essential role in quantum systems. On the one hand, it causes decoherence of quantum states; thus, it limits the coherent dynamics. Therefore, protection of the quantum states from the coupling to the environment has been a crucial issue in quantum engineering. On the other hand, the dissipation can be used as an efficient tool for the preparation and manipulation of particular quantum states of interest (*1*, *2*). Understanding and controlling nonequilibrium dynamics of correlated quantum many-body systems with dissipation are an imperative issue shared in common among experimental systems in diverse areas of physics, including ultracold gases (*1*–*3*), Bose-Einstein condensates (BECs) placed in optical cavities (*3*, *4*), trapped ions (*5*, *6*), exciton-polariton BEC (*3*, *7*), and microcavity arrays coupled with superconducting qubits (*8*, *9*).

Cold atoms, which attract much attention for the investigation of quantum many-body systems owing to high controllability of various parameters, are often regarded as an ideal closed (or isolated) quantum system. However, this controllability also allows the creation of open quantum systems by introducing dissipation processes. So far, various kinds of theoretical works on the effect of dissipation have predicted novel quantum states engineered by dissipation due to photon scattering and particle loss (*10*–*20*). Experimentally, a one-body dissipation has been introduced in a controlled manner with several methods. The utility of an electron beam has been demonstrated in previous studies (*21*–*23*). With a well-designed photon scattering process, measurement backaction on the many-body state (*24*) and the many-body localization in open quantum systems (*25*) have been investigated. In the case of the three-body loss process, controlling the strength of three-body recombination by Feshbach resonance and realization of a novel metastable many-body state have been demonstrated (*26*).

Because the two-body interaction is fundamental and crucial for the emergence of the novel quantum states and many-body physics such as quantum phase transitions, it is important to investigate the effect of two-body dissipation on quantum many-body systems in a manner that the strength of dissipation can be widely controlled. In the reported pioneering works, the two-body loss process was realized by using intrinsic nature of molecules, such as vibrational quenching (*27*) and chemical reaction (*28*), and the lifetime of the molecules was investigated, however a systematic study of the effect of two-body dissipation on quantum many-body physics has not been reported.

Here, we report an investigation of a Bose-Hubbard system using ultracold atoms in a three-dimensional (3D) optical lattice, in which we introduce engineered dissipation of the two-body particle losses. By exploiting the highly controllable nature of the dissipation that we introduce, we successfully reveal the effect of the dissipation on the quantum phase transition from a Mott insulator to a superfluid state in a systematic manner. In particular, we observe in the ramp-down dynamics across the crossover from the Mott insulator to the superfluid states that the melting of the Mott state is delayed and the growth of the phase coherence is suppressed for the strong dissipation. Note that the type of the dissipation introduced in this work is an on-site one. The highly controllable on-site dissipation allows us to study the quench dynamics as the novel method of the initial state preparation, providing a new way for investigating nonequilibrium quantum dynamics. The success in engineering the controllable dissipation of the Bose-Hubbard system offers new opportunities for exploring novel roles of the dissipation in quantum many-body systems.

## RESULTS

### Engineered two-body dissipation

A crucial part of the research of a driven-dissipative quantum many-body state is the design of engineered dissipation for a quantum many-body system. In our experiment, two-body inelastic atom loss with controllable strength is successfully implemented by introducing a single-photon photoassociation (PA) process for ultracold ytterbium (^{174}Yb) atoms in a 3D optical lattice (see Materials and Methods). The PA beam drives the intercombination transition of ^{1}*S*_{0}↔^{3}*P*_{1}, by which two atoms in the doubly occupied sites are photoassociated into the ^{1}*S*_{0} + ^{3}*P*_{1} molecular state and immediately dissociated into the two ground-state atoms (see Fig. 1A and Materials and Methods). This process gives high kinetic energy to the dissociated atoms, which thus results in the escape from the lattice. In this manner, the PA laser induces the two-body inelastic collision loss between the two atoms occupying the same site. We can realize controllable strength of inelastic collision coefficient β_{PA} by changing the intensity of the PA laser *I*, as shown in Fig. 1B. β_{PA} is determined through the relation Γ_{PA} = β_{PA}∫|*w*(** r**)|

^{4}

*d*

**, where the inelastic collision rate Γ**

*r*_{PA}is measured from the exponential fit of collisional loss dynamics, as shown in the inset of Fig. 1B, and

*w*(

**) is the Wannier function of the lowest band. Note that the measurement is done in a deep lattice potential of**

*r**V*

_{0}= 14

*E*

_{R}so that atom tunneling is suppressed in the time scale of this measurement. Here, is a recoil energy, where

*m*is the mass of the

^{174}Yb atom,

*h*is the Planck’s constant, and λ

_{L}= 532 nm is the wavelength of the lattice laser. In this manner, we can realize controllable strength of inelastic collision up to β

_{PA}~ 1.2 × 10

^{−10}cm

^{3}s

^{−1}corresponding to Γ

_{PA}~ 70 kHz in the lattice depth of

*V*

_{0}= 14

*E*

_{R}and γ ~ 5, where γ = ℏΓ

_{PA}/

*U*is the dimensionless dissipation strength that is independent of the lattice depth, and

*U*denotes the on-site interaction.

### Model

It is important to understand how this PA process is effectively described in a dissipative Bose-Hubbard model. Theoretically, the system of bosonic atoms in a sufficiently deep optical lattice coupled coherently to the ^{1}*S*_{0} + ^{3}*P*_{1} molecular state via the PA laser is well described by the Markovian master equation for the coupled atom-molecule mixture model (*29*) with a one-body molecular loss term. By adiabatically eliminating the molecular degrees of freedom based on a second-order perturbation theory for the master equation (*30*, *31*), we derive the effective master equation (see section S1 for details)(1)where(2)and(3)(4) denotes the annihilation operator of atoms at site *j* and . 〈*j*, *k*〉 represents nearest-neighboring pairs of lattice sites. This model is nothing but the single-component Bose-Hubbard model with a two-body loss term (*31*), where *J* and Γ_{PA} denote the hopping energy and the strength of the two-body inelastic collision induced by PA, respectively. In Eq. 4, *g* and Γ_{M} denote the strengths of the atom-molecule coupling and the one-body molecular loss, respectively. Whereas *g* is controllable by varying the intensity of the PA laser, Γ_{M} is fixed for a specific molecular state. Note that the effective master equation (Eq. 1) is valid only when , where *D* and *W* denote the detuning of the PA coupling and the on-site interaction between an atom and a molecule, respectively. Because max(|*g*|, |*D*|, |*W*|, |*U*|, |*J*|)/ℏ ~ 100 kHz at most and Γ_{M} ~ 1 MHz in our experiments (see section S2), this condition is safely satisfied.

### Stability of the atoms with a unit-filling initial state

Before studying the effect of the dissipation on the quantum phase transition, we investigate the stability of the atoms with a unit-filling initial state at a fixed lattice depth. Here, the strength of the dissipation is varied in a wide range from the weak region, in which the dissipation acts as perturbation, to the strong region, in which it exceeds any other energy scale. In contrast to the previous works in which the experiments were done only in the limited range of the dissipation strength, the wide range of our engineered dissipation enables us to observe a crossover between qualitatively different roles of the dissipation. This gives a clue to understand the many-body physics with the dissipation.

In contrast to the measurement of β_{PA} shown in Fig. 1B, this measurement is done at a shallow lattice depth of *V*_{0} = 8 *E*_{R}, in which the tunneling rate 6*J*/ℏ is 4.7 kHz. Thus, in the absence of dissipation, atoms can tunnel to neighboring sites in a time scale of this measurement so that we can investigate how the atom tunneling, which is the only mechanism causing doubly occupied sites, is modified by the dissipation. In our experiment, we first adiabatically load a BEC of 1.0 × 10^{4} atoms into a 3D optical lattice with *V*_{0} = 15 *E*_{R}, in which the state is a singly occupied Mott insulator. Subsequently, we ramp down the lattice to *V*_{0} = 8 *E*_{R} in 0.2 ms and apply the PA laser at the same time. The initial atom number at each site is at most unity, confirmed by the absence of the atom loss by the PA laser and also the occupancy-sensitive high-resolution laser spectroscopy (*32*).

The measured two-body loss rate κ is shown in Fig. 1C. For weak dissipation, the two-body loss rate κ grows as Γ_{PA} increases, which reflects the increase of the detection rate of tunneling. For strong dissipation, however, the two-body loss rate decreases when Γ_{PA} increases. Namely, atom loss is suppressed by the strong on-site dissipation. This counterintuitive behavior is a manifestation of the continuous quantum Zeno effect (*33*), that is, the strong two-body inelastic collision plays a role of the strong measurement and suppresses the coherent process of tunneling. From the comparison between the theory and experiment in a wide range of dissipation strength, we confirm that the measured loss behavior correctly captures the theoretical prediction (see section S3). Note that we observe unexpectedly large atom loss for much higher intensity of PA laser, which prevents the suppression of the two-body loss rate for the strong dissipation region over γ ~ 5 from clear observation (see section S4). Therefore, in our experiment, we restrict the region of the dissipation strength under γ ~ 5.

### Effect of the dissipation on the quantum phase transition

We next investigate the effect of the on-site dissipation on the quantum phase transition from the Mott insulator to the superfluid (*34*), which is the main topic of the present work. Specifically, starting with a singly occupied Mott insulating state, we analyze the dynamics of the atoms subjected to the PA laser during a ramp-down of the lattice depth. The ramp-down speed is −2 *E*_{R}/ms, which is much slower than the case of the two-body loss measurement discussed above. Before presenting the experimental observation, we theoretically analyze such dynamics by assuming a homogeneous system and solving the effective master equation (Eq. 1) within the Gutzwiller mean-field approximation (*15*, *16*) to obtain some insights into the problem. Details of the theoretical analyses are shown in section S3.

An important effect of the dissipation on the quantum phase transition is that it explicitly breaks the conservation of the particle number of the system. Because the superfluid–Mott insulator transition at γ = 0 is originated from the U(1) symmetry associated with the particle number conservation, the introduction of finite γ changes the transition to a crossover. Note, however, that the two-body loss term does not explicitly break the U(1) symmetry. The master equation (Eq. 1) is invariant under the U(1) transformation, , where ϕ is an arbitrary constant.

This crossover can be theoretically characterized by the growth rate of the superfluid order parameter amplitude; when the growth rate is smaller, the system is deeper in the Mott insulator region. According to this characterization, we find that, in the strong on-site dissipation region, where γ ≫ 1, the Mott insulating state is more favored for larger γ (see fig. S10A). This effect originates from the quantum Zeno suppression of the tunneling, which is observed in the two-body loss rate measurement.

This interesting effect of the on-site dissipation on the crossover manifests as the delay in the melting of the singly occupied Mott insulator in the ramp-down dynamics. In Fig. 2A, we show the atom number per site and the condensate fraction as functions of the instantaneous lattice depth during the ramp-down dynamics. We see that, in the strong dissipation region, the onset of the atom loss or the order parameter growth shifts to the side of small lattice depth as γ increases. This result suggests that one may experimentally observe the delay in the Mott insulator melting by measuring the time evolution of the atom number and the momentum distribution during the ramp-down dynamics.

Having the above theoretical insights in mind, we perform the experiment for measuring ramp-down dynamics across the crossover from the Mott insulator to the superfluid. The atom number and the momentum distribution during ramp-down dynamics are obtained from the fluorescence detection and the density distribution of the time-of-flight (TOF) absorption image, respectively. Our experiment starts with ramping up the lattice to *V*_{0} = 20 *E*_{R} for the preparation of the singly occupied Mott insulator state. The atom number is tuned to be small enough that no doubly occupied site exists. Subsequently, we ramp down the lattice by applying the PA laser. The lattice ramp-down speed is −2 *E*_{R}/ms. After ramping down the lattice to the final lattice depth, we performed the fluorescence detection to measure the atom number, or we suddenly turned off all the trap and took the absorption image after 8-ms ballistic expansion to obtain the density distribution.

We first focus on the atom loss measurement during ramp-down dynamics. Figure 2B shows the atom number measured with various dissipation strengths. The experimental result well reproduces the overall features of the calculation shown in Fig. 2A (left). Specifically, the significant atom loss starts around *V*_{0} = 10 *E*_{R} in the presence of weak dissipation (②), whereas the atom number is conserved during ramping down the lattice without dissipation (①). This onset shifts to the deep lattice side as γ increases (③) for weak dissipation (γ < 2). However, when γ increases further from γ ~ 2, the onset shifts to the shallow lattice side (④). To identify the onset, we fit the double linear function to the data (Fig. 2C), which are shown in Fig. 2B. In the presence of on-site dissipation, the atom loss is correlated with the melting of the Mott insulator in ramp-down dynamics because the melting creates the double occupation, which is blasted out by the PA laser. Our result suggests that the melting of the Mott insulator is delayed for strong on-site dissipation. Quantitatively, the onset changes from *V*_{0} = 11.7(4) *E*_{R} to *V*_{0} = 9.2(4) *E*_{R} at the maximum as γ increases. This corresponds to the increase of *zJ*/*U* by a factor of 2.2. These behaviors capture the essence of the theoretical predictions mentioned above.

Figure 3A shows a series of TOF absorption images obtained by changing the final lattice depth from Mott insulator regime to superfluid regime with various strengths of the dissipation. Without dissipation, we observe the transition from a Mott insulator state to a superfluid state, as shown in Fig. 3A of γ = 0: In the deep lattice such as *V*_{0} = 20 *E*_{R}, we obtain a broad distribution with no pattern, which indicates that the atoms have no phase coherence corresponding to the Mott insulator state. By ramping down the lattice across the critical depth of *V*_{0} = 11.3 *E*_{R}, which is calculated from the scattering length of ^{1}*S*_{0} state of ^{174}Yb (*35*), we obtain a clear interference pattern characterizing the presence of the phase coherence of the superfluid state. In the presence of dissipation, the observed transition is significantly modified, as shown in Fig. 3A. As the strength of the dissipation increases, the interference pattern becomes unclear in the shallow lattice regime. For strong dissipation such as γ ~ 5, any pattern cannot be observed. This result indicates that the growth of the phase coherence is suppressed by the strong dissipation.

To quantitatively evaluate the phase coherence, we introduce the visibility of the interference peaks as *v* = (*N*_{max} − *N*_{min})/(*N*_{max} + *N*_{min}) (*36*). Here, *N*_{max} is the sum of the number of atoms in the regions of first-order interference peaks, and *N*_{min} is the sum of the number of atoms in the regions at the same distance from the central peak along the diagonals. Whereas the visibility increases with the ramp-down of the lattice, this increase becomes more moderate in the stronger dissipation, as shown in Fig. 3B. Especially, clear dependence on the strength of the dissipation is observed below the depth of *V*_{0} = 11 *E*_{R}, which is around the calculated critical depth at γ = 0. As shown in Fig. 3B, the effect of the dissipation on the width of the crossover region is observed as more moderate growing of the visibility below the depth of 11 *E*_{R} associated with the increase of γ. In addition, as shown in Fig. 3C, the dissipation moderates the narrowing of the width (in micrometers) of the density distribution, and the slope of the width (in micrometers) with respect to the lattice depth becomes less steep as γ increases. Narrowing the width of the distribution indicates the localization of the state in the momentum space. All these measurements support the delay in the melting of the singly occupied Mott insulator in the ramp-down dynamics as an effect of the on-site dissipation as we see in the calculation of the condensate fraction shown in Fig. 2. Note that the observation of the excitation gap, which is the direct evidence of the formation of the Mott insulator, is difficult in this dissipative system because the excitation spectrum should have a broad linewidth determined by the inelastic collision rate Γ_{PA} of a few tens of kilohertz.

### Quenching the dissipation

It is important to experimentally check whether this behavior is attributed to some heating effect by the PA laser. For this purpose, we measure the phase coherence after turning off the PA laser. If the absence of the interference pattern is attributed to the heating, the phase coherence is no longer restored after the PA laser is turned off. In contrast, if the state after the ramp-down of the lattice is still a Mott insulator, the phase coherence can be restored. Similar to the measurement of Fig. 3A, we ramp down the lattice to a final lattice depth in −2 *E*_{R}/ms with the maximum strength of dissipation γ = 4.6(4). Then, we suddenly turn off the PA laser and investigate the subsequent time evolution of the atoms in the lattice by observing the phase coherence through a TOF absorption image at some hold time.

The result for the case of the final lattice of *V*_{0} = 8 *E*_{R} is shown in Fig. 4B for the observed TOF images and in Fig. 4C for the evolution of the visibility and width of the density distribution. After some hold time, an interference pattern grows. It serves as a direct signature of the restoration of the phase coherence, indicating that the absence of the interference pattern in Fig. 3A is not completely attributed to the heating. We confirm that the total atom number is conserved in this dynamic, as shown in Fig. 4E. This means that the evaporative cooling, which could possibly explain the observed behavior, does not occur during the dynamics. The atoms, after turning off the dissipation, can be considered as an isolated (closed) system. Therefore, the observed dynamics in our experiment should not be considered as the usual thermal relaxation with the environment, but the relaxation in the isolated quantum system, which is a hot topic actively studied in recent experiments and theories (*37*–*40*). Here, we consider the tunneling time as a relevant time scale because the superfluid state is realized through the process of delocalization of the particles by the tunneling. As shown in Fig. 4C, the time constant of the increase of the visibility and the decrease of the width is comparable to the tunneling time (6*J*/*ℏ*)^{− 1} = 0.21 ms. Figure 4D shows the visibility and the width of the density distribution, 0 and 4 ms after the ramp-down for various final lattice depths, similarly indicating the restoration of the coherence.

## DISCUSSION

We have realized the engineered dissipative Bose-Hubbard system by introducing a controllable strength of two-body inelastic collision using a PA laser. By exploiting the highly controllable nature of the dissipation, we have investigated the effect of the dissipation on the quantum phase transition from the Mott insulator state to the superfluid state in the lattice ramp-down dynamics. We have observed that the melting of the Mott state is delayed and the growth of the phase coherence is suppressed for the strong on-site dissipation. The favored state depends on the type of the dissipation. For example, the stabilization of the superfluid state using a well-designed off-site dissipation is proposed (*10*). Because of the marked change in the onset of the Mott insulator melting, as shown in the increase of *zJ*/*U*, we can access the interesting problem of quenching the dissipation across the crossover from the Mott insulator to the superfluid, where turning off the dissipation corresponds to a sudden parameter change of the Bose-Hubbard system (*41*). In this method, the required time for turning off the dissipation could be very short, whereas the sudden change of the depth of optical lattice needs a certain time to stabilize the power of the lattice laser as well as to prevent nonadiabatic interband transition. Moreover, although we have used ^{174}Yb, which is a bosonic isotope of an alkali-earth–like species, to demonstrate our method for controlling the dissipation, it is generally applicable to other atomic species that can be coupled to a state of lossy PA molecule. The crossover properties can also be caused by varying the on-site interaction (see fig. S10B). Controlling the on-site interaction with Feshbach resonance, for example, using alkali atoms, enables us to investigate a wider range of strength of dissipation including infinitely strong, because weakening the on-site interaction corresponds to strengthening the dissipation γ. Our work opens a new way to study the quantum many-body system by controlling the dissipation.

## MATERIALS AND METHODS

### Preparation of ultracold ^{174}Yb atoms in an optical lattice

After collecting atoms with a magneto-optical trap (MOT) with the intercombination transition of ^{1}*S*_{0}↔^{3}*P*_{1}, we loaded the atoms into a crossed far-off resonant trap (FORT). Subsequently, an evaporative cooling was performed, resulting in an almost pure BEC with no discernible thermal component. The atom number of BEC was controlled by changing the collection time of atoms in MOT. The trap frequencies of the FORT at the final stage of the evaporative cooling are = (162, 31,166) Hz. Here, the *x* ′ and *y* ′ axes were tilted from the lattice axes (*x* and *y*) by 45° in the same plane. Then, the BEC was loaded into a 3D optical lattice with the wavelength of the lattice laser λ_{L} of 532 nm. A typical atom number of BEC is 6 × 10^{4} in a measurement of the PA rate of doubly occupied sites, or 1 × 10^{4} in a measurement of two-body loss rate and observation of the quantum phase transition for atoms without multiply occupied sites.

### Details of PA

By applying the PA laser to the atoms in the lattice, two atoms in the doubly occupied sites were photoassociated into the ^{1}*S*_{0} + ^{3}*P*_{1} molecular state whose vibrational quantum number = 16. These molecules were immediately dissociated into the two ground-state atoms. We note that the photon scattering of the atoms due to the PA laser can be negligible because the selected PA line is 3.7 GHz below the atomic transition line, which is about 2 × 10^{4} times larger than its natural line width. β_{PA} for various intensities was determined through the loss dynamics of the atoms by measuring the remaining atom number with the fluorescence detection method after applying the PA laser to the atoms in the lattice with depth of *V*_{0} = 14 *E*_{R}. In this lattice depth, the system is in a state with singly and doubly occupied sites. After the PA laser was applied, the remaining atom number *N*(*t*) decreased as *N*(*t*) = *N*_{1} + *N*_{2} exp(−Γ_{PA}*t*), where *N*_{1} and *N*_{2} are the initial atom number in the singly and doubly occupied sites, respectively.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/12/e1701513/DC1

section S1. Derivation of the dissipative Bose-Hubbard model

section S2. Loss dynamics from the Mott insulating state with double filling

section S3. Details of the theoretical analyses using the Gutzwiller variational approach

section S4. Unexpectedly large atom loss for strong intensity of PA laser

fig. S1. Time evolution of the normalized atom density for = 2.

fig. S2. Measurement of the one-body molecular loss .

fig. S3. On-site interaction *U* and the hopping energy *J* as a function of the lattice depth.

fig. S4. Time sequence of the atom loss measurement from the Mott insulating state with unit filling.

fig. S5. Time evolution of the atom density for .

fig. S6. Loss rate κ as a function of the dissipation strength γ.

fig. S7. Time evolution of ρ_{3,3}.

fig. S8. and as a function of γ.

fig. S9. Time evolution of the amplitude of the superfluid order parameter and its growth rate.

fig. S10. Contour plot of the growth rate, .

fig. S11. Time sequence for the dynamical melting of the Mott insulating state with unit filling.

fig. S12. Atom density as a function of the instantaneous value of the lattice depth *V*_{0}/*E*_{R}.

fig. S13. Condensate fraction as a function of the instantaneous value of the lattice depth *V*_{0}/*E*_{R}.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank Y. Ashida, M. Barbier, and P. Naidon for fruitful discussions.

**Funding:**This work was supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT)/Japan Society for the Promotion of Science (JSPS) KAKENHI grant numbers JP25220711, JP26247064, JP16H00990, JP16H01053, and JP16H00801; Core Research for Evolutional Science and Technology; Japan Science and Technology Agency (JST) no. JPMJCR1673; and Matsuo Foundation. T.T. acknowledges support from the JSPS (KAKENHI grant number JP16J01590).

**Author contributions:**T.T., S.N., and Y. Takasu carried out experiments and the data analysis. I.D. carried out the theoretical calculation. Y. Takahashi conducted the whole experiment. All authors contributed to the writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from T.T. (tomita{at}scphys.kyoto-u.ac.jp).

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).