The problem is discretization error.

density divides the range of data by default into n=512 intervals. The outlier causes the intervals to be at least an order of magnitude greater than otherwise. With very large amounts of data (it seems to take somewhere more than $10^5$ for this problem to become really apparent with the default value of n and this order-of-magnitude outlier), the cumulative discretization error takes over.

The general lesson--which is by no means confined to just density estimation nor to R--is that when routine tools are pushed towards their limits, we need to understand how they do their calculations and we should remember to pay attention to all the parameters that control their accuracy.

I haven't reverse-engineered the source code to pin this down definitively, but I will supply some empirical evidence using random data sets.

trial <- function(n, mean=0, sd=5, outlier=700, kernel="epanechnikov", width=512) {
  set.seed(17) # Makes everything reproducible
  x <- c(outlier, rnorm(n, mean, sd))
  #
  # Plot the densities for visual comparison.
  #
  plot(dd1 <- density(x, kernel=kernel, n=width), col="Red", xlim=c(-sd,sd)*5)
  lines(dd2 <- density(x[-1], kernel=kernel, n=width), col="Blue", xlim=c(-sd,sd)*5)
  #
  # Print crude integral estimates (trapezoidal rule) to memorialize the results.
  # 
  print (c(sum(diff(dd1$x)*(dd1$y[-1] + dd1$y[-length(dd1$y)])/2),
         sum(diff(dd2$x)*(dd2$y[-1] + dd2$y[-length(dd2$y)])/2)))
}
trial(10^6)                       # Looks like the error in the question
trial(10^6, kernel="rectangular") # Has similar results--problem is not the kernel
trial(10^6, width=2^13)           # Both integrals are essentially 1 now
trial(10^3)                       # Agreement is ok
trial(10^3, width=2^6)            # The problem is reproduced

The problem is discretization error.

density divides the range of data by default into n=512 intervals. The outlier causes the intervals to be at least an order of magnitude greater than otherwise. With very large amounts of data (it seems to take somewhere more than $10^5$ for this problem to become really apparent with the default value of n and this order-of-magnitude outlier), the cumulative discretization error takes over.

The general lesson--which is by no means confined to just density estimation nor to R--is that when routine tools are pushed towards their limits, we need to understand how they do their calculations and we should remember to pay attention to all the parameters that control their accuracy.

I haven't reverse-engineered the source code to pin this down definitively, but I will supply some empirical evidence using random data sets.

trial <- function(n, mean=0, sd=5, outlier=700, kernel="epanechnikov", width=512) {
  set.seed(17) # Makes everything reproducible
  x <- c(outlier, rnorm(n, mean, sd))
  #
  # Plot the densities for visual comparison.
  #
  plot(dd1 <- density(x, kernel=kernel, n=width), col="Red", xlim=c(-sd,sd)*5)
  lines(dd2 <- density(x[-1], kernel=kernel, n=width), col="Blue", xlim=c(-sd,sd)*5)
  #
  # Print crude integral estimates (trapezoidal rule) to memorialize the results.
  # 
  print (c(outlier=sum(diff(dd1$x)*(dd1$y[-1] + dd1$y[-length(dd1$y)])/2),
         non.outlier=sum(diff(dd2$x)*(dd2$y[-1] + dd2$y[-length(dd2$y)])/2)))
}
trial(10^6)                       # Looks like the error in the question
trial(10^6, kernel="rectangular") # Has similar results--problem is not the kernel
trial(10^6, width=2^13)           # Both integrals are essentially 1 now
trial(10^3)                       # Agreement is ok
trial(10^3, width=2^6)            # The problem is reproduced
trial(10^3, outlier=10^4)         # The problem is enormous due to the wide range