Christer Hublin
email: christer.hublin@ttl.fi

Markku Partinen

Karoliina Koskenvuo

Karri Silventoinen

Markku Koskenvuo

Jaakko Kaprio

Á M Partinen

Á K Silventoinen

Á M Koskenvuo

Á J Kaprio

Shift-work and cardiovascular disease: a population-based 22-year follow-up study

Keywords: Shift-work, Coronary heart disease, Cardiovascular disease, Mortality, Disability retirement, Genetic effect

come   L'archive ouverte pluridisciplinaire

Introduction

Many work-related psychosocial factors, including work schedules, have been associated with the risk of cardiovascular disease (CVD) [START_REF] Knutsson | Health disorders of shift workers[END_REF]. Despite evidence put forth for a modest excess risk for CVD among shift-workers compared with daytime workers [START_REF] Bøggild | Shift work, risk factors and cardiovascular disease[END_REF] the heterogeneous nature of the studies (17 up to the year 2000) meant that no valid meta-analytic risk estimate could be computed [START_REF] Bøggild | Shift work and heart disease: meta analysis of the epidemiological literature[END_REF]. In a very recent systematic review Frost et al. [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF] concluded that the epidemiological evidence for a causal relation between shift-work and ischemic heart disease (or CHD) is limited. The review is based on 14 studies published before end of March 2008. Seven analyzed fatal events, six combined fatal and non-fatal events, and one study reported separately both types of events. Relative risks ranged from 0.6 to 1.4 in 12 papers while two papers reported relative risks around 2.0. Most studies based on fatal events showed no or weak associations but studies that combined fatal and non-fatal events showed modest positive associations. In a majority of studies negative or positive bias due to the quality of outcome or exposure information, or confounder control could not reasonably be ruled out [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF]. Interestingly, there are only two studies [START_REF] Kawachi | Prospective study of shift work and risk of coronary heart disease in women[END_REF][START_REF] Knutsson | Shiftwork and myocardial infarction: a case-control study[END_REF] assessing the risk in women.

Several hypotheses on the possible causal mechanisms between shift work and CVD or CHD have been presented. These include mismatch of circadian rhythm (affecting nutrition and causing stress related to sleep deprivation, leading to metabolic syndrome), social stress (due to isolation), unfavourable behavioral factors (increased smoking, sedentariness, poor diet), and lower socioeconomic status [START_REF] Bøggild | Shift work, risk factors and cardiovascular disease[END_REF][START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF].

There are few studies on the association between shiftwork and retirement. In a prospective Danish study with 15 years follow-up shift-work was not associated with disability pension in men, but the risk was increased in women even after adjustments for several possible confounders {hazard ratio (HR) = 1.34, 95% confidence interval (CI) 1.02, 1.75 [START_REF] Tu ¨chsen | A 15-year prospective study of shift work and disability pension[END_REF]}. Several Japanese studies have shown an approximately 1.2-fold risk in up to 30% increases both in systolic and diastolic blood pressure [START_REF] Suwazono | Shift work is a risk factor for increased blood pressure in Japanese men: a 14-year historical cohort study[END_REF] and significantly increased risks (1.1-3.6 fold) [START_REF] Morikawa | Relationship between shift work and onset of hypertension in a cohort of manual workers[END_REF][START_REF] Sakata | The relationship between shift work and the onset of hypertension in male Japanese workers[END_REF][START_REF] Oishi | A longitudinal study on the relationship between shift work and the progression of hypertension in male Japanese workers[END_REF] for hypertension have also been reported. Large prospective studies of hypertension risk in shift-workers are lacking in European and US populations.

In this study we assessed the association between type of working time (especially shift-work) in both genders and the risk of mortality due to CHD, disability retirement due to CVD, and incident hypertension in a large, prospective, population-based cohort of Finnish adults with 22-year follow-up data and detailed information on potential confounding or effect-modifying factors. Differences in shift-work history among twin pairs discordant for cardiovascular outcome would support a causal hypothesis between this common exposure and these major public health problems.

Methods

The Finnish Twin Cohort

The Older Finnish Twin Cohort consists of all Finnish twin pairs of the same gender born before 1958 with both co-twins alive in 1975. These pairs were selected from the Central Population Registry of Finland in 1974 [START_REF] Kaprio | Genetic and environmental factors in complex diseases: the Older Finnish Twin Cohort[END_REF]. The first questionnaire was mailed in 1975 (response rate 89%) and the second in 1981 (response rate 84%). The questionnaires included questions on socio-demographics, health status, lifestyle and psychosocial factors including work patterns. The present analyses include those responding to the both questionnaires with information on working time, and resident in Finland in 1981 (N = 20,142; 51.2% women; mean age in 1981 40.2 years).

The study has been approved by the ethical committee of the Department of Public Health, University of Helsinki. Informed consent was obtained from all respondents.

Questionnaire data

Information on working time was asked (both in 1975 and 1981): ''Is your present work or the work you last were doing: (1) mainly daytime work, (2) mainly nighttime work, (3) mainly shift-work, or (4) I have never been working''. Physical work load was enquired about by asking (in 1981): ''What type is/was your work: (1) mainly sitting, (2) standing and walking without other physical activity, (3) in addition to standing and walking also lifting and carrying, or (4) heavy physical work''. Working pace was asked (in 1981): ''How much can you influence your working pace: (1) I can relatively freely regulate it, [START_REF] Bøggild | Shift work, risk factors and cardiovascular disease[END_REF]. I can to some extent regulate it, (3) I can not regulate it, or (4) I am not working''.

We assessed the following socio-demographic and lifestyle covariates (asked both in 1975 and 1981; see Table 1): married or cohabiting, social class (upper or lower white collar, skilled or unskilled workers, farmers, others), education (nine categories by years of school, high school equals to 12 years), smoking status (never, occasional, ex-, or current smoker), binge drinking [START_REF] Kaprio | Genetic influences on use and abuse of alcohol: a study of 5638 adult Finnish twin brothers[END_REF][START_REF] ¨t | Binge drinking in midlife and dementia risk[END_REF], grams of alcohol consumed daily (based on self-reported average quantities of beer, wine and spirits consumed) [START_REF] Romanov | Self-reported alcohol use: a longitudinal study of 12, 994 adults[END_REF], self-report of hypertension (physician-made diagnosis or medication) [START_REF] Hernelahti | Stability and change of volume and intensity of physical activity as predictors of hypertension[END_REF], BMI (body mass index, kg/m 2 ), conditioning physical activity (sedentary, intermediate, vigorous physical activity) [START_REF] Kujala | Relationship of leisure-time physical activity and mortality[END_REF], life satisfaction [START_REF] Koivumaa-Honkanen | Selfreported life satisfaction and recovery from depression in a 1-year prospective study[END_REF], diurnal type [START_REF] Koskenvuo | Heritability of diurnal types: a nationwide study of 9212 adult twin pairs[END_REF], insufficient sleep [START_REF] Hublin | Insufficient sleep: a population based study in adults[END_REF], and sleep length and use of hypnotics and/or tranquillizers [START_REF] Hublin | Sleep and mortality: a population-based 22-year follow-up study[END_REF]. Life satisfaction correlates highly (r = -0.63) with depression as assessed concurrently by the Beck Depression Inventory [START_REF] Koivumaa-Honkanen | Life satisfaction and depression in a 15-year follow-up of healthy adults[END_REF]. Additionally, we used information of working time (mainly daytime, nighttime, two-shift without night shift, two-shift with night shift, three-shift, never been working) from the 1990 questionnaire in a subgroup (N = 10,957, 54.4%) to assess the stability of working time type.

Follow-up data

Vital status (alive in Finland on December 31, 2003, date of death or date of migration from Finland) was obtained from the Population Register Centre of Finland. The follow-up for mortality was from the exact date of response (date 1981 questionnaire returned) to December 31, 2003. Cause-of-death statistics of the underlying cause [due to CHD; ICD-10 codes I20-I25, and corresponding sections of ICD-8 and ICD-9 (codes 410-414)] were obtained from Statistics Finland. Information on disability retirement due to CVD [ICD-10 codes I00-I99, and corresponding sections of ICD-8 and ICD-9 (codes 390-459)] was obtained from the Finnish Social Insurance Institution (FSII) and the Finnish Centre for Pensions. Information of incident hypertension was obtained by record linkage of the Cohort data with the register of FSII on new cases of reimbursement of cost for medical treatment of hypertension. Reimbursement is granted after submission of a medical certificate by the treating physician, reviewed by the medical experts of the FSII. These registers cover all Finnish citizens. Our analyses are based on 402,266 person-years, and 857 deaths due to CHD, 721 disability retirements caused by CVD, and 2,642 new cases of hypertension during the follow-up.

Data analysis and statistical methods

The study population was grouped according to the type of working time (each subject belonging to one group) as follows (Table 1 shift-work both in 1975 and 1981. These categories were used as our main interest was to compare the risks associated with constant shift work (group 5) to that with constant daytime work (group 1). Those ever working night-time formed a separate group as their number was small and they differed in many characteristics from the other groups (Table 1). Those not working or with missing data on working status were excluded (N = 1,096) from the analyses. Subjects on disability pension or retired prior to 1981 were excluded from the analyses of risk of disability retirement, and similarly those with baseline hypertension from the analysis of incident hypertension during the followup 1982-2003. Significance of differences in the prevalence of each covariate between the types of working time was assessed by computing likelihood-ratio tests from logistic regressions models including age with and without the shiftwork categories. As there are well known gender differences in CHD and CVD, the results are given separately for men and for women.

Cox proportional hazard models were used to obtain HR and their 95% CI for mortality due to CHD and disability retirement caused by CVD by type of working time. As the reference group we used those daytime working both in 1975 and 1981. The proportional-hazard assumptions were tested and confirmed using Schoenfeld residuals by examining that the log hazard-ratio function is constant over time for the shift-work variables, using the phtest option in Stata, which was used in all statistical analyses (version 9.2, Stata Corporation, College Station, Texas, USA). Also, plots of the -ln(-ln(survival) vs. ln(age) were used to visually inspect the proportional hazards assumption. Because the study sample included twin pairs, standard errors and CIs were adjusted for possible within-pair correlations using robust estimators of variance [START_REF] Williams | A note on robust variance estimation for clustercorrelated data[END_REF].

Additionally, as the follow-up time is long, we also computed HRs for the mortality risk due to CHD up to 1990 to exclude the possibility that the risk might accumulate non-uniformly.

The risk of mortality and disability retirement was also calculated for the healthy-subgroup of the total study sample (N = 15,314). This includes the subjects without chronic diseases such as angina pectoris, cancer, myocardial infarction, stroke, diabetes mellitus or chronic obstructive pulmonary disease in 1981 based on the questionnaire and information from multiple medical registers [START_REF] Kujala | Modifiable risk factors as predictors of all-cause mortality: the roles of genetics and childhood environment[END_REF][START_REF] Sørensen | Intention to lose weight, weight changes, and 18-y mortality in overweight individuals without co-morbidities[END_REF]. The subjects with hypertension only were not excluded.

In the fully adjusted models, adjustments were made for the socio-demographic and lifestyle covariates (measured in 1981) that are potential confounders associated both with shift-work and risk of death (see above). The subjects with missing data on any of the covariates were excluded from the fully adjusted models.

Additionally, our twin data permitted us to compute polychoric correlations, which assess pairwise similarity of those shift-working both in 1975 and 1981 and those daytime working both in 1975 and 1981. The analysis was based on 916 male and 1,081 female monozygotic (MZ) pairs, and 1950 male and 2006 female dizygotic (DZ) pairs. Discordant pair analysis permits assessment of the possible causal relationship between shift-work and disease. This was performed in those pairs with one twin dying due to CHD during the follow-up (218 MZ and 568 DZ pairs), to see if these twin-pairs differed in shift-work pattern while adjusting for latent familial factors shared by the twins, including partially (DZ pairs) or completely for genes (MZ).

Results

Demographic characteristics are given in Table 1. Both in 1975 and 1981, 76.9% were daytime working, 9.5% shiftworking, and 1.7% had night-time work at either point. In 1975, the most common job categories among men were metalworker, worker in process industry and working in transportation, and in women nursing, factory worker and assistant jobs in nursing/health care. In the subgroup with information of the working time type in 1990, of those with daytime work both in 1975 and 1981, 90.8% were still daytime working, and of those shift-working both in 1975 and 1981, 68.8% were still shift-working (of them 59.4% had night shifts).

When comparing the characteristics between different types of working time, there were significant differences among both genders in marital status, social class, education, work load and pace, smoking, diurnal type, short sleep, and insufficient sleep. Significant differences were found in physical activity and in life satisfaction among males, and in BMI and in long sleep among females. Between daytime working and shift-working, there were some notable differences. The proportions of low social class and low level of education were higher among shiftworkers. Among daytime workers there were less of those unable to regulate working pace and more of those with high physical work load. There were more short sleepers and evening-type persons among shift-workers, but these proportions were even higher among night-time workers. The proportion of those with higher income level was highest among shift-workers. Low life satisfaction was most prevalent among night-time workers and least prevalent among permanent shift-workers. The proportion of those frequently using hypnotics and/or tranquillizers was on the same level in all groups.

Initiating smoking between 1975 and 1981 (i.e. became current smokers) was significantly associated with working-time types only in women. Those women who had started shift-work after 1975 had a 1.5-fold greater risk (HR = 1.49, 95% CI 1.03, 2.16) of smoking initiation compared to consistent daytime workers. However, when former smokers were excluded from the analysis, the association attenuated and became non-significant, suggesting that initiating shift-work may be associated with a relapse to smoking among former smoking women.

The associations of type of working time with the risk of CHD mortality and with the risk of disability retirement due to CVD are shown in Table 2. The age-adjusted HRs for CHD mortality were highest but non-significant in night-time workers. No significant increase in the risk of CHD mortality was found among permanent shift-workers. In the full model the HRs were similar in the total population (with the exception of night-time working women), and in the healthy-subpopulation the risk was even lower. Expectedly an elevated risk of CHD mortality was found for current smokers, for those with hypertension and for those men with obesity (BMI [ 25 kg/m 2 ). Thus, there was no significantly increased risk of CHD mortality among shift-and night workers as compared to permanent daytime workers. Also the estimates for mortality up to 1990 were essentially similar.

Age-adjusted risk of disability retirement due to CVD was significantly increased only in women starting shiftwork (HR = 1.88; 95% CI 1.10, 3.20). In the full model, both in the total population and in the healthy-subpopulation, the HRs were generally decreased, with the statistically non-significant exceptions of healthy men stopping shift-work and healthy women starting shift-work. Low social class, hypertension and obesity (BMI [ 25 kg/m 2 ) were associated with a higher risk of disability retirement due to CVD. In addition, current smoking and use of hypnotics/tranquillizers increased the risk of disability retirement in men (data not shown in Tables).

There were no significant increase in the risk of incident hypertension either in age-adjusted or in the fully adjusted model.

We carried out subgroup analyses. The CHD mortality risk by age-groups was analyzed separately for the both genders showing several increased HRs in a non-consistent pattern. No significantly increased risks were found except for middle-aged men who had started shift-work in 1981 (HR = 2.09, 95% CI 1.14-3.86). After full adjustment, the HR decreased and became non-significant.

In pairs discordant for CHD, the odds ratios for shiftwork were non-significant, indicating no difference in shift-working patterns between twins in pairs where one consequently died from CHD.

Discussion

The present results do not support a direct association between shift-work and cardiovascular (CVD), especially coronary heart disease (CHD) risk, neither in men or women. There was no increase in the risk of mortality due to CHD, disability retirement due to CVD, or incident hypertension among shift-workers compared to daytime workers during the 22-year follow-up. The results were similar also in the healthy subpopulation free of chronic disease. In twin pairs discordant for CHD, there was also no evidence for a shift-work effect.

Our study has several strengths: the large study sample (totally about 400,000 person-years) is representative of the general population with very high response rates, the number of cases was substantial (857 deaths and 721 retirements), the type of working time was measured twice, the information on potential confounding or effect-modifying factors was detailed, and the statistics on causes of deaths and retirement are comprehensive and reliable.

Additionally, in our study information also on women was provided. There are no indications that rates of CHD and blood pressure levels would differ between singletons and twins [START_REF] Christensen | Cardiovascular mortality in twins and the fetal origins hypothesis[END_REF][START_REF] De Geus | Comparing blood pressure of twins and their singleton siblings: being a twin does not affect adult blood pressure[END_REF], and there is no reason why the association between shift work and outcome should be different in adult twin individuals and singletons doing similar work. The recent systematic review [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF] assessed 14 studies reporting on the risk of ischemic heart disease (or CHD) in relation to shift-work. Eight of these reported mortality data. In six studies the relative risks were 0.64-1.11, based on 163-8,767 cases. An increased risk was found only in one study [START_REF] Kawachi | Prospective study of shift work and risk of coronary heart disease in women[END_REF] with female American nurses (relative risk 1.19, 95% CI 0.63-2.23, N = 79,109) and in another study [START_REF] Fujino | Study Group. A prospective cohort study of shift work and risk of ischemic heart disease in Japanese male workers[END_REF] with male Japanese workers (HR = 2.32, 95% CI 1.37, 3.95, N = 17,649) but the number of cases were low (44 and 86, respectively). Two studies used employed males from general populations [START_REF] Fujino | Study Group. A prospective cohort study of shift work and risk of ischemic heart disease in Japanese male workers[END_REF][START_REF] Virtanen | Socioeconomic inequalities in cardiovascular mortality and the role of work: a register study of Finnish men[END_REF] and five studies included male workers from industrial settings. The followup period was [15 years in four studies, and confounder control was limited in three studies [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF]. In two studies [START_REF] Kawachi | Prospective study of shift work and risk of coronary heart disease in women[END_REF][START_REF] Fujino | Study Group. A prospective cohort study of shift work and risk of ischemic heart disease in Japanese male workers[END_REF] the information on working hours was based on self-report. In seven studies the source of outcome has been based on death certificates. Methodologically, the present study is comparable, in some respects superior to the previously published mortality studies.

There is a striking difference in our results compared to those of Fujino et al. [START_REF] Fujino | Study Group. A prospective cohort study of shift work and risk of ischemic heart disease in Japanese male workers[END_REF], despite methodological similarities: large population-based prospective study with long follow-up, self-report of working hours, adequate adjustment for possible confounders, and similar amount of person-years (233,869 vs. in the present study 193,061 male person-years). However, there may be variation in the CHD/CVD patterns between populations, explaining at least partly this difference. Also the socioeconomic status of shift-workers and exposure to other risk work-related and other factors may differ between countries: traditionally in Finland a typical male shift-worker is a well-earning paper-mill worker, whereas a typical female shift-worker works in low-income jobs in healthcare and social work sectors. Especially, there was a large difference in the number of CHD deaths (in Fujino et al. 86 CHD-deaths in the whole cohort and of them 18 in shift-workers, in our study 536 and 53, respectively). Serum levels of cholesterol and the prevalence and incidence of CHD are considerably lower in Japan compared to Western counterparts [START_REF] Ueshima | Explanation for the Japanese paradox: prevention of increase in coronary heart disease and reduction in stroke[END_REF], which might affect the results.

Hypertension, dyslipidemia and smoking are the major acknowledged risk factors of CHD and CVD. However, it has been suggested that adjustment for these is not appropriate as they should not be considered confounders but rather steps in the causal chain [START_REF] Bøggild | Shift work, risk factors and cardiovascular disease[END_REF][START_REF] Tu ¨chsen | A 12 year prospective study of circulatory disease among Danish shift workers[END_REF][START_REF] Knutsson | Shift work and ischaemic heart disease[END_REF]. Doing so would make it difficult to compare results in shift-work with other studies. There is evidence that smoking is more prevalent among shift-workers than daytime workers [START_REF] Bøggild | Shift work, risk factors and cardiovascular disease[END_REF][START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF][START_REF] Van Amelsvoort | Smoking among shift workers: more than a confounding factor[END_REF]. Our results suggest that initiating shift-work may be associated with a relapse to smoking among former smokers (especially in women), which is in accordance with the fact that smoking is generally initiated in adolescence before entering the work-life. As the relationships are complex, giving both adjusted and unadjusted results is advisable [START_REF] Knutsson | Methodological aspects of shift-work research[END_REF]. In the present study adjustment did not essentially change the results, not supporting a risk of over-adjustment. We found no association of shift-work with hypertension, suggesting it cannot be a step in the putative chain from shift-work to cardiovascular outcomes, while adjustment for smoking did not affect the estimates related to shift-work. Lipid measures were not available in our study.

Regarding other outcomes our results differed from the findings of previous studies. Overall, we found no significantly increased risk of disability pensions in women as reported in the Danish study [START_REF] Tu ¨chsen | A 15-year prospective study of shift work and disability pension[END_REF], although age-adjusted risk was increased in women starting shift-work. We did not either found increased risk of hypertension as in Japanese studies [START_REF] Morikawa | Relationship between shift work and onset of hypertension in a cohort of manual workers[END_REF][START_REF] Sakata | The relationship between shift work and the onset of hypertension in male Japanese workers[END_REF][START_REF] Oishi | A longitudinal study on the relationship between shift work and the progression of hypertension in male Japanese workers[END_REF].

Some limitations in our study need to be considered. First, the question in 1975 and 1981 on the type of working time did not specify the presence of night shifts. More detailed information of the details of the shift system would have been desirable as different systems may not be associated with similar effects, and a broad definition of shift-work could attenuate the association to CVD. Methodologically, a detailed assessment of the working schedules would be challenging in a study population derived from general population, as there is a multitude of different shift systems. However, studies assessing effects of shift work with and without night work, have not found clear indications for stronger effects among shift workers mainly working during nights compared to those who did not work at night [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF]. In Finland in 1984, 44% of men and 35% of women shiftworking were working in three shifts [START_REF] Lehto | Threats and Opportunities[END_REF]. It is also unclear what relation is between the risk and the amount of years of shift-work [START_REF] Frost | Shift work and the risk of ischemic heart disease-a systematic review of the epidemiologic evidence[END_REF]. In the present study shift-working seemed to be long-lasting as in a subgroup with information of working time in 1990 about 2/3 were still working in shifts of those who were shift-working both in 1975 and 1981. Secondly, it may be argued that disability retirement due to CVD is not a reliable outcome as it is possible that at least a part of the workers with CVD change to daytime work rather than become disability retired, which might dilute our result. However, this is not supported by a Finnish study among nurses showing that the possibility of employees with cardiovascular problems leaving their workplace was independent of work schedule [START_REF] Kivima ¨ki | Prevalent cardiovascular disease, risk factors and selection out of shift work[END_REF]. Third, we did not study a cohort of shift workers but a sample of general population with subjects with different types of working hours. Therefore, this study is free from risks which may be associated with a single line of business, and our results may be applicable to the population in general. However, both approaches have their relevant aspects.

In conclusion, our results do not support an association between shift-work and increased risk of CVD, CHD, or hypertension either in men or in women. More prospective studies with large representative populations, long followups, reliable measure of exposure to shift-work, and adequate control of confounders are needed.

Table 1

 1 Selected demographic characteristics in 1981 of the entire cohort with working time data in 1975 and 1981

		P-value a (m/f)			NA			\0.01/0.02	\0.01/\0.01	\0.01/\0.01	\0.01/0.01	\0.01/0.01	\0.01/\0.01	\0.01/\0.01	0.71/0.32	0.12/0.04	0.73/0.74	0.05/0.85	\0.01/0.08	\0.01/\0.01	\0.01/\0.01	0.16/\0.01	\0.01/\0.01	0.94/0.58
		Shift-work both	in 1975 and 1981	N (m/f) = 962/930	495/505	339/280	128/145	78.7/64.3	78.0/59.0	2.7/8.6	66.3/19.4	5.4/4.3	18.0/17.2	39.9/21.7	16.5/5.4	43.6/22.6	13.2/12.7	18.3/7.3	13.9/14.5	11.8/11.6	18.9/17.2	12.5/18.1	18.8/23.6	2.1/2.5
		Daytime work	in 1975, shift-work in 1981	N (m/f) = 562/610	426/443	100/115	36/52	65.0/66.0	71.6/55.3	9.3/14.5	41.6/7.4	7.4/5.0	13.8/15.4	43.6/29.5	15.6/7.4	32.8/22.2	10.3/11.6	17.1/7.7	18.0/14.3	13.5/12.3	23.3/14.5	12.4/18.3	19.8/23.5	1.9/3.3
		Shift-work in 1975,	daytime work in 1981	N (m/f) = 509/701	349/520	132/115	28/66	66.6/68.0	68.0/53.0	10.8/19.2	50.9/10.3	8.8/3.3	8.3/8.7	39.9/24.0	17.2/5.0	35.8/18.8	11.0/10.6	21.0/8.8	20.8/17.1	10.1/10.6	18.6/13.3	15.4/25.9	16.0/19.7	2.6/2.7
		Night-time either	in 1975 or 1981	N (m/f) = 108/224	64/136	28/60	16/28	69.9/70.1	56.4/65.4	13.3/8.2	43.9/14.0	4.5/4.8	8.9/11.6	46.0/31.2	20.4/6.0	39.8/26.4	15.7/11.6	13.9/7.1	22.1/18.5	24.3/23.4	30.6/24.5	13.5/16.7	24.3/33.8	2.4/2.0
	Type of working time	Daytime both in	1975 and 1981	N (m/f) = 7,698/7,838	4,159/4,270	2,296/2,078	1,243/1,490	73.9/67.5	52.6/41.2	13.1/16.1	53.3/16.8	15.6/6.8	3.8/7.5	35.6/19.4	15.2/6.0	39.5/25.2	12.5/14.5	17.1/7.0	16.0/16.1	8.7/8.6	14.0/12.7	14.6/20.6	13.4/21.2	2.7/3.0
				Age	24-39 years	40-54 years	55 years or more	Married	Social class: bluecollar workers	Education: high school or more	Higher income level b	High physical load of work	Unable to regulate working pace	Current smoker	Heavy alcohol use c	BMI C 25	Hypertension	Conditioning physical activity	Low life satisfaction	Diurnal type: evening	Short sleeper (B7 h)	Long sleeper (C8 h)	Insufficient sleep d	Frequent use of hypnotics/tranquillizers e

Percentage in each working time type of those in the specified category of potential confounder; number of men/women

(9,

839/10,303, total N = 20,142) given in column header. Likelihoodratio tests age-adjusted NA not applicable a Likelihood-ratio test of difference between working time groups b Percentage of those with yearly income exceeding the lower limit (1800 Finnish Marks in 1975) of the best earning one-third of the study population (genders combined) c On average, more than a drink a day for women and more than 2 for men d A difference C1 h between the self-reports of sleep need and sleep length e C60 days/year of either medication Shift-work and cardiovascular disease

  ): group 1. daytime work both in 1975 and 1981, group 2. night-time work either in 1975 or 1981, group 3. shift-work in 1975 and daytime work in 1981, group 4. daytime work in 1975 and shift-work in 1981, and group 5.

Table 2

 2 Risk of mortality due to coronary heart disease (CHD), risk of disability retirement due to cardiovascular disease (CVD) and risk of incident hypertension in 1982-2003 [hazard ratioIn addition to age includes as covariates marital status, social class, education, smoking status, binge drinking, grams of alcohol consumed daily, hypertension, body mass index (BMI), conditioning physical activity, life satisfaction, diurnal type, sleep length, use of hypnotics and/or tranquillizers, physical load of work, and working pace

			Shift-work both	in 1975 and 1981			1.09 (0.82, 1.44)	1.22 (0.83, 1.79)		1.06 (0.75, 1.50)	1.21 (0.75, 1.93)	0.77 (0.42, 1.39)	1.02 (0.27, 3.95)			1.15 (0.84, 1.58)	0.96 (0.60, 1.55)		0.70 (0.48, 1.03)	0.79 (0.43, 1.43)	0.88 (0.56, 1.37)	0.82 (0.40, 1.70)			1.15 (0.97, 1.37)	1.06 (0.88, 1.28)		1.07 (0.88, 1.30)	1.00 (0.80, 1.23)
	(95% confidence interval)]	Type of working time	Daytime both Night-time work either Shift-work in 1975, Daytime work in 1975,	in 1975 and 1981 in 1975 or 1981 daytime work in 1981 shift work in 1981	Mortality due to CHD	Age-adjusted	Men (N = 9,839, n = 536) 1.00 1.75 (0.92, 3.33) 1.09 (0.68, 1.76) 0.94 (0.56, 1.56)	Women (N = 10,303, n = 321) 1.00 1.38 (0.71, 2.69) 1.16 (0.66, 2.04) 1.55 (0.88, 2.74)	Full model a	All men (N = 9,175, n = 458) 1.00 1.82 (0.97, 3.41) 0.86 (0.48, 1.54) 0.79 (0.44, 1.41)	All women (N = 9,434, n = 250) 1.00 0.90 (0.36, 2.23) 1.08 (0.56, 2.08) 1.52 (0.82, 2.82)	Healthy men (N = 7,146, n = 152) 1.00 1.17 (0.35, 3.90) 0.14 (0.02, 1.05) 0.91 (0.39, 2.12)	Healthy women (N = 6,936, n = 36) 1.00 NC b 0.67 (0.09, 4.96) 2.16 (0.62, 7.55)	Disability retirement due to CVD	Age-adjusted	Men (N = 8,298, n = 403) 1.00 0.89 (0.34, 2.32) 1.30 (0.82, 2.04) 1.46 (0.90, 2.38)	Women (N = 8,719, n = 219) 1.00 1.68 (0.82, 3.46) 1.03 (0.55, 1.95) 1.88 (1.10, 3.20)	Full model a	All men (N = 7,831, n = 375) 1.00 0.51 (0.17, 1.49) 1.12 (0.70, 1.78) 1.13 (0.67, 1.90)	All women (N = 8,073, n = 188) 1.00 0.94 (0.35, 2.51) 0.93 (0.48, 1.82) 1.32 (0.72, 2.40)	Healthy men (N = 6,711, n = 245) 1.00 0.63 (0.15, 2.60) 1.43 (0.83, 2.44) 1.32 (0.70, 2.49)	Healthy women (N = 6,565, N = 126) 1.00 1.03 (0.32, 3.27) 0.80 (0.32, 2.00) 1.84 (0.94, 3.60)	Incident hypertension	Age-adjusted	Men (N = 9,474, n = 1,305) 1.00 0.73 (0.39, 1.35) 0.93 (0.72, 1.21) 1.04 (0.81, 1.33)	Women (N = 9,775, n = 1,337) 1.00 1.12 (0.80, 1.56) 1.05 (0.84, 1.32) 1.24 (0.99, 1.56)	Full model a	All men (N = 8,841, n = 1,242) 1.00 0.72 (0.37, 1.39) 0.93 (0.71, 1.21) 1.04 (0.81, 1.33)	All women (N = 9,000, n = 1,215) 1.00 1.02 (0.70, 1.48) 0.97 (0.76, 1.23) 1.15 (0.89, 1.47)	On each row given number of subjects (N) and number of cases (n) included in the analysis. Statistically significant associations given in bold

a b Not computable

C. Hublin et al.

Acknowledgments Supported by the Academy of Finland Center of Excellence in Complex Disease Genetics and the Social Insurance Institution of Finland.


Felix Schmitt

Markus Hild

Robert Roth

Ab initio phase diagram of ultracold 87 Rb in a one-dimensional two-color superlattice

Keywords: numbers: 67.85.Hj, 03.75.Lm, 67.85.-d

We investigate the ab-initio phase diagram of ultracold 87 Rb atoms in an onedimensional two-color superlattice. Using single-particle band structure calculations we map the experimental setup onto the parameters of the Bose-Hubbard model. This ab-initio ansatz allows us to express the phase diagrams in terms of the experimental control parameters, i.e., the intensities of the lasers that form the optical superlattice. In order to solve the manybody problem for experimental system sizes we adopt the density-matrix renormalizationgroup algorithm. A detailed study of convergence and finite-size e ffects for all observables is presented. Our results show that all relevant quantum phases, i.e., superfluid, Mott-insulator, and quasi Bose-glass, can be accessed through intensity variation of the lasers alone. However, it turns out that the phase diagram is strongly affected by the longitudinal trapping potential.

Introduction

Ultracold atomic gases in optical lattices have been a topic of active research for about a decade now. One of the research thrusts is the use of these systems as experimental quantum simulators for a variety of lattice models and allow for detailed investigations of strongly correlated quantum systems in a perfectly controllable environment [1,2]. By tuning the laser intensities of the optical lattice alone, one can seamlessly drive a system through quantum phase transitions like the superfluid to Mott-insulator transition [3,4]. In so-called two-color superlattices, additional lasers are used to introduce irregular lattice topologies which give rise to exotic quantum phases like a quasi Bose-glass phase [5,6].

Strongly correlated particles in periodic potentials are well described by Hubbard-type models. Together with powerful many-body methods this allows for theoretical studies of the phase diagram of ultracold atomic gases in optical lattices [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. However, a one-to-one comparison between experiment and theory has rarely been done so far because these theoretical studies usually adopt the generic parameters of the Hubbard model to span the phase diagram. Such a phase diagram of ultracold bosonic atoms in a two-color superlattice is shown in Fig. 1(a).

In this work we establish a closer link to experiments by computing the phase diagrams with respect to the natural experimental control parameters, which are the intensities ( s 2 , s 1 ) of the two lasers generating the one-dimensional optical superlattice. Such an experimentspecific phase diagram is shown in Fig. 1(b). In order to predict this type of phase diagram we start with single-particle band structure calculations to extract the Hubbard parameters for a specific experimental setup. Then, the many-body problem is solved using the densitymatrix renormalization group (DMRG) algorithm. Together with quantum Monte Carlo (QMC) methods [25], the DMRG is one of the most powerful many-body techniques for the treatment of one-dimensional lattice problems. After introducing and benchmarking our band structure plus DMRG approach, we discuss the phase diagram in the ( s 2 , s 1 )-plane for different The labels mark the domains of the superfluid (SF) phase, the homogeneous Mott-insulator (MI) phase, and the quasi Bose-glass (BG) phase (taken from Ref. [22]).

observables in this strictly one-dimensional framework. The laser wavelengths, longitudinal and transverse trapping frequencies, and atomic species are chosen from a range covering recent experimental setups [4,5,6]. We focus on the dependence of the phase diagrams on the transverse trapping frequency ω ⊥ and the longitudinal trapping frequency ω x , which was not discussed in detail before.

1D Bose-Hubbard Model and Band Structure Calculations

The single-band Bose-Hubbard model [26] is a widely used framework for studying the ultralow temperature physics of strongly correlated, neutral atoms in su fficiently deep optical lattices. We assume a one-dimensional lattice with I sites and N bosonic atoms. For each site we define the creation (annihilation) operators â † i (â i ) with respect to the localized Wannier states corresponding to the lowest Bloch band. The mean occupation-number at each lattice site is given by ni = â † i âi . The Bose-Hubbard Hamiltonian

Ĥ = I i=1 -J i,i+1 ( â † i+1 âi + â † i âi+1 ) + 1 2 U i (n i -1) ni + i ni (1)
accounts for three basic processes: the tunneling of atoms to adjacent sites, the on-site twobody interaction, and the on-site potential energy. The site-dependent Hubbard parameters J i,i+1 , U i , and i define the relative strengths of the individual terms and contain all information about depth and topology of the optical potential, and the interaction between the atoms. The phase diagrams spanned directly by these parameters, typically using J i,i+1 ≡ J, U i ≡ U and some ansatz for i to account for superlattice structures, are extensively discussed in Refs. [12,13,14,15,16,19,20,21]. More recently, the on-site energies i were calculated directly from the parameters of the optical superlattice to provide a closer connection to experiment [23].

In this work, our aim is a discussion of the phase diagram using the experimental parameters directly and not the generic Hubbard parameters. To this end, an explicit treatment of the underlying single-particle physics is necessary. Therefore, we start from the optical potential generated by two orthogonal polarized standing-wave laser-fields with wavelengths λ 1 and λ 2 and the respective potential depths s 1 and s 2 . Furthermore, we consider an additional harmonic potential with frequency ω x accounting for the intensity variation of the optical lattice through the focusing of the laser beams and a magnetic trapping potential. Using the recoil energy E r i = h 2 2mλ 2 i of atoms with mass m as a natural energy scale and a phase shift φ between the standing waves, the potential along the x-axes reads:

V(x) = s 1 E r 1 sin 2 2π λ 1 x + φ + s 2 E r 2 sin 2 2π λ 2 x + 1 2 mω 2 x x 2 . ( 2
)
Throughout this work we consider a setup defined by λ 2 = 800 nm and s 2 for the primary laser generating the optical lattice potential and λ 1 = 1000 nm and s 1 for the secondary laser generating the two-color superlattice topology with a phase shift of φ = π/4.This leads to the commensurate superlattice that was also used in previous publications [14,15,19,20,21,22].

Before we are able to extract the Hubbard parameters for a given potential, we have to determine the localized Wannier functions via a single-particle band structure calculation. For a periodic potential (s 1 = 0, ω x = 0) with I sites, we numerically obtain the solutions for the Bloch functions ψ k (x) in the lowest energy band. The quasimomenta are quantized with respect to the size I of the optical lattice and are labeled by k with k = 0, 1, • • • , I -1. A Fourier transformation of the Bloch functions with respect to the quasimomenta in the subspace of the lowest energy band leads to the Wannier functions

w i (x) = 1 √ I I-1 k=0 ψ k (x)e -i 2π I ki e iϕ k . (3)
The arbitrary phases ϕ k are chosen such that the resulting Wannier functions are maximally localized at their individual lattice site i. Using these maximally localized Wannier functions, the Hubbard parameters are obtained via the matrix elements of the individual terms of the real-space Hamiltonian [27]

-J i, j = dx w * i (x) - 2 2m ∂ 2 ∂x 2 + V(x) w j (x) i = dx w * i (x) - 2 2m ∂ 2 ∂x 2 + V(x) w i (x) ( 4 )
U i = 2 ω ⊥ a s dx |w i (x)| 4 .
The contact interaction is defined by the three-dimensional s-wave scattering length a s . The transverse directions are integrated out assuming Gaussian wavefunctions with frequencies ω y = ω z = ω ⊥ [2]. This one-dimensional description is valid as long as the tunneling in the transverse directions is strongly suppressed, i.e., as long as the laser intensities in these directions are sufficiently large. Following Refs. [5,6] we consider a gas of 87 Rb atoms with s-wave scattering length a s = 109 r Bohr and we assume a transverse trapping frequency ω ⊥ = 2π × 17 kHz. In the last part of this work we will discuss the changes in the phase diagram induced by a different value of ω ⊥ . As a first application of our band structure calculations we validate the single-band approximation in the Hubbard model. As a by-product from the calculations of the Hubbard parameters we obtain from the single-particle band structure calculation the energy gap ∆E between the first and the second band of Bloch functions. In Table 1 we list some values of ∆E together with the Hubbard parameters U and J in the relevant parameter range. Since ∆E is always about one order of magnitude larger than U and J, we conclude that excitations to energetically higher Bloch bands induced by tunneling or interaction can be neglected even for shallow optical lattices.

As a second application of our band structure calculations we check for the validity of the restriction to nearest-neighbor tunneling and on-site two-body interactions. For the weakly and the strongly interacting regime, we calculate the respective matrix elements of the Hubbard Hamiltonian using Eqs. (4) and U i, j = 2 ω ⊥ a s dx |w i (x)| 2 |w j (x)| 2 for the interaction term. The results are shown in Table 2. Even in the weakly interacting regime (s 2 = 2), the nearest-neighbor tunneling exceeds more-distant tunneling processes by at least one order of magnitude. The interaction matrix element for neighboring lattice sites is already two orders of magnitude smaller than the on-site interaction matrix element. In the strongly interacting regime (s 2 = 10) we already have two orders of magnitude between J i,i+1 and J i,i+2 and five orders of magnitude between U i,i and U i,i+1 . Since we focus on the intermediate and strong interaction regime, the restriction to J i,i+1 and U i is well justified. So far we have discussed the limit of a homogeneous optical lattice. As soon as the secondary laser which generates the superlattice, or an additional harmonic potential are taken into account, a straight-forward band structure calculation is no longer doable, because Bloch functions are only defined for strictly periodic potentials. Therefore, in order to extract site-dependent Hubbard parameters also for an inhomogeneous lattice we are limited to an approximate scheme to obtain localized Wannier functions. We use two di fferent approaches to extract the site-dependent Hubbard parameters.

As a simple ansatz, we consider the secondary laser as a perturbation of the strong primary laser (s 1 s 2 ). The Wannier functions are extracted from a conventional band structure calculation for a homogeneous lattice defined by the primary laser alone. In this approximation the Wannier functions are identical for each lattice site. Using these Wannier functions the Hubbard parameters of each site of the superlattice are computed. The site-dependence of the parameters thus results exclusively from the superlattice potential V(x) entering into the matrix elements (4) and not from a site-dependence of the Wannier functions themselves. As a result, the parameter U i characterizing the on-site interaction remains constant for all lattice sites. An exemplary set of site-dependent Hubbard parameters calculated in this scheme is shown in Fig. 2. Please note that we always subtract a global energy constant from the Hamiltonian to set min = min{ i } = 0. In a more sophisticated scheme we determine the site-dependent Wannier functions individually for each site of the inhomogeneous lattice using a standard band structure calculation for a periodic lattice with a lattice amplitude defined by the local depth of the inhomogeneous potential at that particular site. In this way, the shape of the Wannier functions depends nontrivially on the local structure of the superlattice potential. The only reason why the set of Wannier functions determined in this way cannot be considered as an exact set of localized basis functions results from the minimal violation of orthogonality for the Wannier functions of neighboring sites. Their mutual overlap is nonzero but always below 1% in the parameter regime considered in all our calculations. Using these individual localized Wannier functions all site-dependent Hubbard parameters are computed without further approximations. An exemplary set of results is also shown in Fig. 2.

The comparison of the site-dependent Hubbard parameters i and J i,i+1 obtained by the two schemes shows very little difference. This leads to the conclusion that the second and superior scheme provides a sufficiently accurate description of the Hubbard parameters in the parameter range under consideration.

The dominant effect on the Hubbard parameters induced by the superlattice is the spatial variation of the on-site energies i . This is in agreement with the approximation of the superlattice through this parameter alone [14,15,19,20,21]. However, the tunneling matrix element J i,i+1 , which essentially depends on the potential barrier between sites i and i + 1, also shows a significant variation. The on-site interaction matrix element U i exhibits only a weak variation which is introduced by the site-dependence of the Wannier functions in our second scheme.

Density-Matrix Renormalization-Group

We solve the many-body problem associated with the Bose-Hubbard Hamiltonian via the density-matrix renormalization-group (DMRG) algorithm [28,29] which is among the most powerful quasi-exact methods available for one-dimensional lattice models. The so-called infinite-size DMRG algorithm is based on an iterative growing procedure. The algorithm is schematically depicted in Fig. 3. The individual steps are: (0) We start with a block composed of I b sites and up to N b particles described in a Fock space F b of dimension D b . Since the Hubbard Hamiltonian conserves the particle number, the matrix representation of the block Hamiltonian has a block-diagonal form. Each block of the matrix corresponds to a Hilbert space with a fixed particle number. (1) To the block we attach an additional lattice site with up to N s particles described in a Fock space F s to build the Fock space of the system F sys = F b ⊗ F s with dimension D sys = D b D s . Again the matrix representation of the system Hamiltonian is block diagonal. (2+3) In order to simulate a larger lattice, the system is coupled to an analogously constructed environment yielding the superblock H super = F sys ⊗ F env of dimension D super which is projected to a fixed total particle number, satisfying N/I = 1 in our case. All operators coupling the system to the environment (which will later couple the new block to the new site) and all observables have to be transformed accordingly. This cycle is repeated until the final length of the lattice is reached.

The key feature of this algorithm lies in the use of the eigenvectors with largest eigenvalues of the reduced density-matrix as a new, truncated basis for the new block. One can show that this procedure yields an optimized wavefunction, gives the best approximation to expectation values of observables, and preserves a maximum of entanglement between system end environment [29].

The error in the DMRG algorithm is due to the loss of information during the non-unitary basis transformations. It can be estimated by summing up the eigenvalues of the discarded eigenvectors. A smaller sum consequently means a smaller loss of information. In addition one has to consider the restriction to a maximal number of particles max {n i } per lattice site. In the complete Hilbert space this would be equal the total number of particles N. In general, the stronger the correlations between the particles, i.e., the larger Ū, the smaller the sum of the residual eigenvalues and the better the approximation.

If disorder is introduced, then only at the very last step of the growing procedure the full information about the superlattice topology is available to the Hamiltonian. This leads to a poor approximation of the ground state when using the infinite-size algorithm only. As an improvement the finite-size DMRG is applied. After a complete run of the infinite-size algorithm up to the desired length of the lattice, the length of the superblock is kept fixed and the system grows on the expense of the environment and vice versa. During a back and forth sweeping, the superlattice topology is sampled while the Hamiltonian always takes the whole lattice into account. The sweeping continues until all observables are converged.

The way the transformation matrices O are constructed is not uniquely defined by the DMRG algorithm. We would like to emphasize that the reduced density-matrix ρ red is block diagonal and each block has a well defined particle number. One can either use the eigenvectors with the largest eigenvalues for each subspace of ρ red , or one can strictly use the first D b eigenvectors with the largest eigenvalues not accounting for the block-diagonal structure of ρ red . In the first scheme one might discard eigenvectors with sizable eigenvalues if the respective subspace has reached its preassigned dimension. In the second scheme one might discard a complete subspace of a certain particle number in case it has no eigenvector with corresponding eigenvalue among the largest D b eigenvalues. If in a subsequent step of the finite size algorithm this subspace becomes important again, this might prevent the algorithm from converging to the proper ground state. This can be overcome by adding noise to the transformation matrices O during the first few sweeps in the finite-size algorithm with the aim of recovering lost subspaces again [30]. As a third strategy one can keep at least one or a few eigenstates from each subspace even if their eigenvalues are not among the largest D b eigenvalues. We employ the first strategy because it is technically very convenient. However, we checked individual eigenspectra of ρ red and confirmed that none of the discarded eigenvectors had sizable eigenvalues.

Observables

In this section, we introduce the set of observables we employ to distinguish the di fferent quantum phases.

Maximum Number Fluctuation.

The number fluctuation at lattice site i is given by the variance of the occupation number

σ i = n2 i -ni 2 . ( 5
)
The number fluctuation provides information about the local mobility of the atoms in the optical lattice. In order to reduce the amount of information we only consider the maximum number fluctuation through the lattice

σ max = max{σ i } . ( 6
)
Condensate Fraction. In order to determine the fraction of atoms that undergo Bose-Einstein condensation we adopt the Onsager-Penrose criterion [31] and calculate the natural orbitals via the eigensystem of the one-body density-matrix ρ (1) i j = â † i â j . The largest eigenvalue N c of the one-body density-matrix is associated with the number of condensed atoms and defines the condensate fraction

f c = N c N . (7)
Visibility. In the experiment, most information about the atoms in the optical lattice is extracted from the interference pattern obtained by the time-of-flight method. The interference pattern I(δ) is intimately connected to the quasi-momentum structure of the many-body state and can be calculated from the Fourier transformation of the one-body density-matrix [12]

I(δ) = 1 I I i, j=1 e i(i-j)δ ρ (1) i j . (8)
The visibility of the interference fringes ν is obtained from the maxima and minima of the interference pattern

ν = max{I(δ)} -min{I(δ)} max{I(δ)} + min{I(δ)} . ( 9
)
Energy Gap. Measuring the excitation spectrum of the system also provides a sensitive probe for the different quantum phases. In the experiment one employs two-photon Bragg spectroscopy via an intensity modulation of the optical lattice [4,6]. The width of the central interference peak is used as a measure of the energy transfer into the atomic cloud [4]. The detailed structure of the excitation spectrum has been investigated theoretically [19,20,17,18]. Basic information about the excitation spectrum is given by the energy gap ∆E, which is the minimum amount of energy needed to excite the system. It is defined by the difference between the energy of the first excited state and the ground state An inherent complication in the DMRG framework is the calculation of observables. This is because a DMRG calculation does not yield eigenstates of the Hamiltonian in a simple occupation-number basis representation which could be used to compute observables directly. Rather, the matrix representations of all observables have to be dragged through all the cycles of the DRMG algorithm, i.e., they have to undergo all the lossy non-unitary basis rotations. This is why we thoroughly test our DMRG calculations for convergence.

∆E = E 1 -E 0 . (10)

Benchmark of the DMRG algorithm

Convergence

Before we employ the DMRG algorithm to compute phase diagrams for realistic lattice sizes and particle numbers, we have to assess the precision of the numerical DMRG results. We follow a twofold strategy.

First, we compare results for the various observables obtained by DMRG calculations with results from an exact diagonalization scheme [14,15] for a small system with I = N = 10, where the latter calculations are feasible. In the phase diagram shown in Fig. 1(a) we observed an error of the DMRG calculation below 1% for all observables at U/J > 3 already for a small DMRG basis with dimension D super = 338. The complete Hilbert space used in the exact diagonalization scheme has a dimension of D = 92378 for the I = N = 10 system.

Second, in order to validate the results of our DMRG calculations for larger lattices, where no exact calculations in the complete Hilbert space are available, we study the dependence of the DMRG results on basis sizes and particle number truncations used in the algorithm [28,29]. If the results for all observables do not change while the bases size is increased further, the calculation is converged to the exact result. The di fferent basis sets we employ are summarized in Table 3, where max {n i } is the maximum number of particles per lattice site included in the basis. For all following calculations we applied three sweeps in the finite-size algorithm. For simplicity we consider straight lines through the parameter plane shown in Fig. 1(a).

Superfluid to Mott-insulator ( max = 0). The upper row of images in Fig. 4 shows the observables across the superfluid to Mott-insulator phase transition calculated using the three bases specified in Tab. 3. Since the DMRG algorithm is tailored to describe strongly correlated systems, we expect better agreement of the three di fferent calculations with increasing U/J. Apart from the energy gap this is confirmed by our calculations. Only for U/J 3 we observe small differences between the calculations for σ max and f c . For all values of U/J the energy gap is slightly larger when employing the DMRG-A basis. This is because we do not explicitly target at the first excited state for the calculation of the energy gap. Although the ground state has already converged even for the small DMRG-A basis, the first excited state needs a larger basis to converge as well.

Mott-insulator to quasi Bose-glass (U/J = 30). We already pointed out the importance to use the finite-size DMRG algorithm in order to obtain a converged ground state especially when irregularities in the optical lattice are considered. The results of the observables through the Mott-insulator to quasi Bose-glass transition depicted in the lower panel of Fig. 4 show that the finite-size algorithm is perfectly converged for all values of max /J already for the DMRG-A basis.

Since, in this manuscript the focus is on the regime of intermediate and strong interactions, we conclude from our findings that already the DMRG-A basis is suitable to approximate all observables with sufficient precision. Nevertheless, we decided to use the larger DMRG-B basis for all following calculations. Calculations based on this basis are still numerically feasible on a desktop PC while providing good results also in the weakly interacting limit.

Size dependence

We also have to address the dependence of the observables observables on the size of the system. Current experiments typically have between 1 .5 • 10 4 to 2 • 10 5 atoms in the optical lattices [5,6], i.e., roughly between 25-60 atoms in each one-dimensional array of sites. Thus, we will compare DMRG calculations for I = N = 10, I = N = 30, and I = N = 50 all using the DMRG-B basis. In analogy to the above discussion, we show plots through the superfluid to Mott-insulator transition and the Mott-insulator to quasi Bose-glass transition.

Superfluid to Mott-insulator ( max /J = 0). The results are shown in the upper row of Fig. 5. By definition, the maximum number fluctuation is a local observable which is calculated at one individual lattice site and is, therefore, practically independent of the size of the lattice. The energy gap as well as the visibility show only small differences between the small and the two larger lattices indicating a minor dependence on length of the lattice for those observables. However, the condensate fraction depends systematically on the size of the lattice. The larger the lattice is, the steeper is the decrease of f c around U/J ≈ 5. One can easily show that for U/J → ∞ and I = N the condensate fraction scales like f c ∝ 1/I [12] which is in-line with our calculations. We also performed an additional calculation for the large I = N = 50 lattice using the DMRG-C basis. These results are not shown in the plots because there are no sizable deviations to calculations with the DMRG-B basis. Only for the condensate fraction at U/J < 3 the DMRG-C basis yields slightly larger values, e.g. f c = 0.68 instead of f c = 0.63 at U/J = 1. This indicates the slower convergence of the DMRG algorithm in the weakly interacting regime. The results for all other observables remain essentially unchanged when going to the larger DMRG-C basis. Mott-insulator to quasi Bose-glass (U/J = 30). The lower row of Fig. 5 reveals that the energy gap as well as the maximum number fluctuation do not change with the size of the lattice across the Mott-insulator to quasi Bose-glass transition. The condensate fraction exhibits the previously mentioned 1/I scaling which is characteristic for large values of U/J. The visibility in the small lattice is again slightly smaller compared to the two larger lattices.

Considering this analysis, we conclude that calculations including N = 30 particles on I = 30 lattice sites are sufficient to describe realistic experiments. Firstly, because this system size is right in the experimental range. And secondly, when going to larger systems, there are only small and predictable changes for the condensate fraction whereas all other observables remain unchanged.

Ab-Initio Phase Diagrams

After the validation of our framework we now discuss the experiment-specific phase diagram motivated by the experiments [5,6]. We consider an ultracold 87 Rb gas with scattering length a s = 109 r Bohr in an optical lattice with wavelength λ 2 = 800 nm. The superlattice topology is generated by an additional laser with wavelength λ 1 = 1000 nm and relative phase shift of φ = π/4. The respective optical potential depth resulting from the two lasers are given by the dimensionless parameters s 2 and s 1 . The remaining transverse lasers of the optical trap enter via the transverse trapping frequency ω ⊥ which is chosen to be 2π × 17 kHz. Initially, the longitudinal trapping frequency ω x is set to 0 Hz.

We have already used these parameters in Fig. 1 to compare the experiment-specific phase diagram spanned by s 2 and s 1 with a generic phase diagram spanned by U/J and max /J neglecting the site dependence of U and J. Both panels of Fig. 1 show the energy gap ∆E for I = N = 30 obtained from a DMRG calculation using the DMRG-B basis. Since the variation of s 2 and s 1 affects all Hubbard parameters simultaneously, the ( s 2 , s 1 ) phase diagram is distorted in comparison to the (U/J, max /J) phase diagram. However, the (s 2 , s 1 ) phase diagram reveals that all relevant quantum phases are accessible through the variation of the intensity of the two longitudinal lasers alone, while keeping the other parameters fixed.

A detailed analysis of the phase diagram for this set of parameters is given in Figs. 6(a)-(d), where we show the energy gap, the condensate fraction, the maximum number fluctuation, and the visibility, respectively.

The superfluid (SF) phase is characterized by a vanishing energy gap, large condensate fraction, large number fluctuations, and maximum visibility. Although we do not compute the most stringent order parameter for the SF phase -the superfluid fraction [8,12,23] -the above signatures allow us to identify the SF phase in the region of small s 2 up to s 2 6 in the whole range of s 1 shown here. Due to the shallow optical potential in this region the tunneling term in the Hubbard Hamiltonian (1) dominates. This results in a coherent many-body state which is a prerequisite for the SF phase. For s 2 = 6 along 0 < s 1 ≤ 2 the mean interaction energy is Ū/ J ≈ 4.5 which explains the presence of the SF phase in the whole range of s 1 .

In a homogeneous lattice ( s 1 = 0 or max /J = 0) a transition from the SF phase to the homogeneous Mott-insulating (MI) phase occurs around U/J ≈ 5 [23,12] which corresponds to s 2 = 6.25. This is in-line with our results, because around s 2 ≈ 6 the energy gap steeply increases while the condensate fraction, the number fluctuations, and the visibility decrease. At s 2 = 16 and s 1 = 0 the ratio of U/J is 60 and the system is deep in the homogeneous MI phase showing the characteristic large energy gap and vanishing number fluctuations, condensate fraction, and visibility.

If we now increase s 1 at fixed s 2 = 16, the modulation of the site-dependent Hubbard parameters grows rapidly and at s 1 ≈ 0.6 the spread of the on-site energies becomes comparable to the average interaction energy, i.e., max / J ≈ Ū/ J. Thus, despite the strong repulsive interaction, it becomes advantageous to move an atom from a site with large on-site energy to an already occupied site with small on-site energy. Due to this redistribution of particles the homogeneous MI phase is broken up and the transition to the quasi Bose-glass (BG) phase occurs. The commensurate superlattice defined by λ 2 = 800 nm, λ 1 = 1000 nm, and φ = π/4 exhibits only 5 different on-site energies. This small set of on-site energies leads to extended domains in the phase diagram. Two of these domains are visible in Fig. 6(a). Only in the transition region between them the energy gap vanishes.

We emphasize that the quasi Bose-glass phase emerging in a commensurate two-color superlattice is not a gapless phase like the genuine Bose-glass, which only occurs in an infinite lattice with truly random on-site energies. Therefore, we use the term quasi Bose-glass. Intuitively the genuine Bose-glass results from a continuous distribution of on-site energies permitting the construction of excited states by moving particles to sites with infinitesimally larger on-site energies associated with infinitesimally small excitation energies. We have approached this limit using an incommensurate superlattice in a previous publication [22].

Longitudinal Trapping Frequency ω x

Since the aim of this manuscript is the calculation of an experiment specific phase diagram for a realistic experimental setup, it is compulsory to consider an additional magnetic trapping potential and the intensity variation of the optical lattice through the focusing of the laser beams. To this end we have introduced a harmonic potential with frequency ω x in Eq. ( 2). Typical experimental parameters range from ω x = 2π × 10 Hz to 2π × 75 Hz [4,5,6].

To get a impression of the energy scales, we show some values for the Hubbard parameters obtained by our band structure approach in Tab. 4. By setting s 1 = 0 the onsite energies are solely due to the additional harmonic potential. At the outer rims of the lattice (sites 1 and 30) they have the value max .

Up to ω x = 2π × 25 Hz, max / J is an order of magnitude smaller than Ū/ J. For this reason, the phase diagram remains practically unaltered between ω x = 0 Hz and 2π × 25 Hz as can be seen by comparing Figs. 6(a)-(d) and 7(a)-(d). For ω x = 2π × 50 Hz Tab. 4 shows max / J is still about a factor 3 smaller than Ū/ J. As a consequence the onset of the BG phase in Figs. 7(e)-(h) already appears at s 1 ≈ 0.4 instead of s 1 ≈ 0.6 for ω x = 0 Hz. Besides the earlier onset of the BG phase also its gross structure changes. The lobe around s 1 = 0.6 in Fig. 7(e) is suppressed compared to the calculations for ω x < 2π × 50 Hz. Also the maximum fluctuations indicate that the redistribution of particles becomes smoother. This is because for ω x = 0 Hz the superlattice topology exhibits only 5 different on-site energies. With the additional harmonic potential the number of di fferent on-site energies increases and, therefore, the extended domains in the BG phase shrink and become fragmented. For ω x = 2π × 75 Hz the parameters max / J and Ū/ J become comparable and the phase diagram changes dramatically. In Fig. 7(i) the homogeneous MI domain shrinks to a small region (s 2 = 12 -16 and s 1 = 0 -0.2). Furthermore, a clear detection of the BG phase becomes difficult since the characteristic increase of the visibility along the MI to BG transition is no longer visible in Fig. 7(l). Table 4 reveals that at ω x = 2π × 100 Hz the on-site energies max clearly dominate the energy scale. Thus, even the transition from the SF phase to the MI phase is no longer observable in the investigated parameter range of s 2 . The MI domain has completely disappeared in Fig. 7(m) and the visibility remains large throughout the whole range of s 2 for s 1 < 0.2. We note that the uneven contour lines visible in the plots of the maximum number fluctuation σ max for the larger values of ω x have a physical origin. First of all, σ max is particularly sensitive to changes in the intrinsic structure of the many-body state. With increasing ω x the BG phase becomes more fragmented, i.e., as function of the control parameters the number states dominating the ground state change rapidly and a ffect the maximum fluctuations. Moreover, since we plot the maximum of the set of all on-site number fluctuations across the lattice, neighboring points in the phase diagram can show the on-site number fluctuations for completely different lattice sites. Already this will generate kinks in the contour lines. We have checked by increasing the DMRG basis size that this is not an artifact of the DMRG basis truncation.

From the above discussion we conclude that the smaller the longitudinal trapping frequency ω x , the easier is a clear distinction between the SF, MI and BG phases. Thus, any experiment with a focus on the phase diagram of ultracold atoms in an optical superlattice should be designed such that the longitudinal trapping frequency is kept small.

Transverse Trapping Frequency ω ⊥

Finally we study the dependence of the ( s 1 , s 2 ) phase diagram on the intensity of the transverse lasers through a variation of the transverse trapping frequency ω ⊥ . For the sake of simplicity we assume ω x = 0 Hz.

From Eq. ( 4) it follows that the interaction energy U i is proportional to ω ⊥ while i and J i are independent of ω ⊥ . A larger value of ω ⊥ will, therefore, shift the SF to MI transition towards smaller s 2 . Also the spread of the on-site energies must increase to overcome energy cost of a double occupancy and consequently the MI to quasi BG transition will shift towards larger s 1 . In Fig. 8 we show the phase diagrams for ω ⊥ = 2π × 40 kHz and all other parameters unchanged. The gross structure of the phase diagram remains the same. However, in accordance with our considerations, the energy gap is ∆E/ J = 60 already at s 2 = 12.5 instead of s 2 = 16 for ω ⊥ = 2π × 17 kHz (both s 1 = 0). Furthermore, the transition from the homogeneous MI to the quasi BG phase occurs around s 1 ≈ 1.5 as compared to s 2 = 0.6 for ω ⊥ = 2π × 17 kHz.

As a consequence of the dependence of the ( s 2 , s 1 ) phase diagram on the transverse trapping frequency, the precise identification of the phase boundaries is intimately connected to a well defined value of ω ⊥ .

Summary & Conclusions

We have studied the experiment-specific phase diagram of ultracold 87 Rb atoms in an onedimensional two-color superlattice with respect to the parameters of the experiment. Band structure calculations were employed to obtain the generic parameters of the Hubbard model from the experiment-specific parameters. These band structure calculations were also used to confirm the applicability of the Hubbard model in the investigated parameter range.

In order to solve the many-body problem for realistic lattice lengths and particle numbers we have used the density-matrix renormalization-group algorithm. Through a thorough benchmark of our DMRG calculations we demonstrated that all observables are perfectly converged and can practically be considered as exact solutions of the many-body problem. Furthermore, a detailed finite-size analysis for all observables has underlined the significance of our results for realistic experimental system sizes.

Our calculations of the phase diagrams show that all relevant quantum phases can be accessed by only varying the intensities of the two lasers that generate the optical superlattice. For a longitudinal trapping frequency ω x < 25 Hz all different quantum phases can be clearly distinguished by means of the presented observables. However, larger values of the longitudinal trapping frequency lead to radical changes in the structure of the phase diagram and make a clear identification of the quasi Bose-glass phase impossible.

We also showed that the gross structure of the phase diagram does not depend on the transverse trapping frequency ω ⊥ , i.e., the intensity of the lasers in the directions perpendicular to the 1D lattice. However, due to the linear dependence of the interaction energy on the transverse trapping frequency, the position of the transition lines in the phase diagram crucially depend on that parameter.

Figure 1 .

 1 Figure 1. (color online) Contour plots of the energy gap computed using DMRG for a commensurate superlattice with I = N = 30 as (a) a function of the generic Hubbard parameters and (b) a function of the experimental laser intensities.The labels mark the domains of the superfluid (SF) phase, the homogeneous Mott-insulator (MI) phase, and the quasi Bose-glass (BG) phase (taken from Ref.[22]).

Figure 2 .

 2 Figure 2. Site-dependent Hubbard parameters obtained from band structure calculations for a two-color superlattice. Simple ansatz (gray symbols) and calculations with individual Wannier functions (black symbols), both for s 2 = 10 and s 1 = 1 with λ 2 = 800 nm, λ 1 = 1000 nm, ω ⊥ = 2π × 17 kHz, a s = 109 r Bohr , and mass m of 87 Rb. Lines to guide the eye.

Figure 3 .

 3 Figure 3. Sketch of the DMRG cycle. The black dot in step 5 marks the additional site that was attached without increasing the dimension of the Hilbert space of the block. For details see text.

( 3 )

 3 The ground state | ψ 0 is obtained by diagonalizing the superblock Hamiltonian where one can exploit the sparseness of the Hamilton matrix and use efficient Lanzcos or Jacobi-Davidson algorithms. (4) The reduced density-matrix is formed by tracing out the environment ρred = Tr env | ψ 0 ψ 0 | . (4+5) The D b eigenvectors of the reduced density-matrix for the largest eigenvalues are used to span the Fock space for a new block Fb of length Ĩb = I b + 1. These eigenvectors build a non-unitary transformation matrix O which is employed to construct the new block Hamiltonian Hb = O † H sys O.

= 30 Figure 4 .

 304 Figure 4. Benchmark of the convergence of the DMRG calculations. From left to right: energy gap ∆E, maximum number fluctuation σ max , condensate fraction f c , and visibility ν for a I = N = 30 lattice. Upper panels: SF-MI transition at fixed max /J = 0. Lower panels: MI-BG transition at fixed U/J = 30. All plots show three lines corresponding to a di fferent basis choice: DMRG-A (dotted), DMRG-B (dashed), and DMRG-C (solid).

= 30 Figure 5 .

 305 Figure 5. Finite size analysis for I = N = 10 (dotted), I = N = 30 (dashed), and I = N = 50 (solid) based on the DMRG-B basis. From left to right: energy gap ∆E, maximum number fluctuation σ max , condensate fraction f c , and visibility ν. Upper panels: SF-MI transition at fixed max /J = 0. Lower panels: MI-BG transition at fixed U/J = 30.

Figure 6 .

 6 Figure 6. (color online) Phase diagram in terms of energy gap ∆E/ J, condensate fraction f c , maximum number fluctuation σ max , and visibility ν for I = N = 30, a s = 109 r Bohr , ω ⊥ = 2π × 17 kHz, and ω x = 0 Hz.

Figure 7 .

 7 Figure 7. (color online) Phase diagram in terms of energy gap ∆E/ J, condensate fraction f c , maximum number fluctuation σ max , and visibility ν for I = N = 30, a s = 109 r Bohr , ω ⊥ = 2π × 17 kHz. First row: ω x = 2π × 25 Hz, second row: ω x = 2π × 50 Hz, third row: ω x = 2π × 75 Hz, fourth row: ω x = 2π × 100 Hz.

Figure 8 .

 8 Figure 8. (color online) Same parameter set as in Fig. 6 but with transverse trapping frequency ω ⊥ = 2π × 40 kHz. While the gross structure of the phase diagram is independent of ω ⊥ , the scales of the s 2 and s 1 axes change.

Table 1 .

 1 Calculations of the Hubbard parameters U and J and the energy gap ∆E between first and second Bloch band for a homogeneous lattice ( s 1 = 0, ω x = 0) with λ 2 = 800 nm, ω ⊥ = 2π × 17 kHz, a s = 109 r Bohr , and mass m of 87 Rb .

	s 2	2	4	6	8		10	12	14	16
	∆E/E r 2	1.00	1.97	2.90	3.77	4.57	5.31	5.97	6.58
	U/E r 2	0.143 0.205 0.237 0.261 0.280 0.297 0.311 0.323
	J/E r 2	0.143 0.086 0.051 0.031 0.019 0.012 0.008 0.005
	U/J	1.138 2.407 4.671 8.480 14.62 24.22 38.84 60.64
			s 2 = 2 (U i,i /J i,i+1 = 1.38)
	|i -j|	0	1	2		3		4	5
	J i, j /E r 2	-	0.1428	-0.02		0.0048 -0.0014 ≈ 10 -4
	U i, j /E r 2 0.162 0.0035	≈ 10 -4	≈ 10 -5	≈ 10 -6	≈ 10 -7
			s 2 = 10 (U i,i /J i,i+1 = 14.62)
	|i -j|	0	1	2		3		4	5
	J i, j /E r 2	-	0.0192 ≈ -10 -4 ≈ 10 -6 ≈ -10 -8 ≈ 10 -9
	U i, j /E r 2	0.28 ≈ 10 -5	≈ 10 -8	< 10 -12 < 10 -12 < 10 -12

Table 2 .

 2 Higher order tunneling and interaction energies for a homogeneous lattice ( s 1 = 0, ω x = 0) with λ 2 = 800 nm, ω ⊥ = 2π × 17 kHz, a s = 109 r Bohr , and mass m of 87 Rb .

Table 3 .

 3 Bases used for studying the convergence of the DRMG calculations. See text for details.

	D b	D s	D super max{n i }
	DMRG-A 21	126	446	5
	DMRG-B 56	336	5073	5
	DMRG-C 210 1470 68356	6

Table 4 .

 4 Comparison of the Hubbard parameters to analyze e ffect of the longitudinal trapping potential. The parameters are: λ 2 = 800 nm, λ 1 = 1000 nm, φ = π/4, s 1 = 0, ω ⊥ = 2π × 17 kHz, and mass and scattering length of 87 Rb.

	s 2	2	10	12	16
	J/E r 2	0.1428	0.0192	0.0123	0.0053
	Ū/E r 2	0.1624	0.2505	0.2966	0.3232
	Ū/ J	1.1378	14.623	24.222	60.636
			ω x = 2π × 25 Hz
	max /E r 2 0.02431 0.02431	0.02431 0.02431
	max / J	0.17052 1.26739	1.98446 4.55600
			ω x = 2π × 50 Hz
	max /E r 2 0.09737 0.09737	0.09737 0.09737
	max / J	0.68196 5.06927 7.937412 18.2389
			ω x = 2π × 75 Hz
	max /E r 2 0.21909 0.21909	0.21909 0.21909
	max / J	1.53407 11.4048	17.8576 41.0344
		ω x = 2π × 100 Hz
	max /E r 2 0.38950 0.38950	0.38950 0.38950
	max / J	2.72634 20.27242 31.7430 72.9423

Acknowledgment

Supported in part by the ExtreMe Matter Institute EMMI in the framework of the Helmholtz Alliance HA216/EMMI and the DFG Sonderforschungsbereich 634.